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NOMENCLATURE 


Field SI 
B - volume factor ft3/scf m3/ms?> 
Cc - compressibility factor psia! kPa-l 
5 - integration constant 
d - inside pipe diameter in cm 
D - inside pipe diameter ft m 
h - reservoir thickness ft m 
k - permeability md um? 
M - molecular weight lbm/Ib-mole g/g-mole 
) - porosity % % 
9 - gas density bm/ft3 kg/m3 
Pp - pressure psia kPa 
q - flow rate STB/D m3/s 
r - radius ft m 
t - time S S 
S - saturation % % 
at - temperature R K 
U - viscosity cp mPa-s 
ny - Boltzmann Transform 
Z - compressibility factor dimensionless 
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Subscript 


es - critical property 

D - dimensionless quantity 
e - external border 

i - initial 

j - number reference 

e) - oul 

r - reduced property 

Ww - well bore 

wa - water 
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ABSTRACT 


The purpose of this study was to examine the effect of the existence of acid 
gases and the use of each of two viscosity correlations (Carr, Kobayashi and 
Burrows and Dean and Stiel) in the determination of the Pseudo Pressure term . 
Three situations were analysed: the first included the use of all available corrections 
for the presence of acid gases; the second ignored these corrections; the third 
involved an analysis based upon gas which was acid free. The same approaches 


were taken using actual gas well test data. 


These situations were examined for a group of gases covering the entire 
pressure range for each of the correlations used. The same study was repeated for 
actual field tests to determine the effects of the various methods and correlations on 


the interpretation of the tests. 


For the theoretical situation, i.e., the full reduced pressure range, parallel 
results of viscosity and Pseudo Pressure were calculated for the various treatments of 
the acid gases. When the results of Pseudo Pressure were calculated from actual tests 
data, the use of the three acid gas treatments and the two viscosity correlations had a 


negligible effect on the interpretation of Pressure Transient Analysis test results. 


* Pseudo Pressure - The Pseudo Pressure term (‘Y¥ = 2f(P/uz) dP) is 
primarily used in the calculation of gas flow behaviour in porous media. 
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CHAPTER I 


INTRODUCTION 


Understanding gas properties, and predicting their behaviour with greater 
accuracy are two of the goals of research in the petroleum field. In addition to the 
economic benefits which result from improvement in the ability to predict reservoir 
size, to determine maximum production rates, and to size field equipment and 
processing plants, these studies may also yield a better understanding of reservoir 
behaviour which will enhance resource management. Pressure Transient Analysis 


(PTA) is an important tool for understanding the reservoir and its behaviour. 


In PTA of gas reservoirs, three methods - Pressure, Pressure Squared, and 
Pseudo Gas Pressure - are used in solving the Darcy-Continuity equation. The first 
two methods, the Pressure and Pressure Squared approaches, require simplifying 
assumptions to solve a non-linear, second order, partial differential equation. Al- 
Hussainy and Ramey (1966, pp 624-642) did not make these assumptions and 


introduced the third approach, the Pseudo Pressure solution. 


The purpose of this study was to investigate the evaluation of the Pseudo 


Pressure term by carrying out the following: 


1. Determine whether or not the assumption of constancy of P/uz and Uz is 


correct and, if so, to determine if there are any limits. 
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Evaluate the cu product required in the evaluation of a compressible 
gas. 

Evaluate the effects of CO2 and H2S on P/uz, 1/1z and Uc terms. 
Examine the differences between the Lee, Gonzalez, and Eakin (1966), 
the Carr, Kobayashi, and Burrows (1954), and the Dean and Stiel 
(1965) viscosity correlations. 

Compare predicted Pseudo Pressures using the Dean and Stiel viscosity 
correlation to those predicted using the Carr et al. viscosity correlation. 
Examine the effects of these viscosity correlations on the interpretation 


of field tests. 
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LITERATURE REVIEW 


1.1 COMPRESSIBILITY FACTOR FOR NATURAL GAS 


An equation of state (EOS) relates a known quantity of gas to fluid density at 
a given or known pressure and temperature; fluid density is inversely proportional to 
volumetric behaviour at a given or known pressure and temperature. For real gases, 
the deviation of volumetric behaviour from its ideal behaviour is expressed by the gas 
compressibility factor (z-factor). There are several approaches available to predict 
this behaviour. These include models based on van der Waals forces such as the one 
proposed by Peng and Robinson (1976), as well as a popular empirical method 
developed by Standing and Katz (1940). 


Although the Standing and Katz method was readily accepted from 1940 
until the mid 1970s for all natural gases, it is not particularly suitable for computer 
applications. Therefore, several mathematical models were created to represent the 
chart developed by Standing and Katz. Three of the more popular models are 
Starling (1975), Hall and Yarborough (1973), and Dranchuk et al. (1974). The 
initial Dranchuk et al. (1974) model, which used only eight variables, was later 
modified in a paper by Dranchuk and Abou-Kassem (1975), to increase the accuracy 
of the z-factor for reduced temperatures below 1.05. Using the modified eleven- 
variable EOS, the Dranchuk and Abou-Kassem (1975) correlations for the z-factor, 


became more similar to the correlations presented by Starling and Hall's EOS. 
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However, the overall performance of the Dranchuk and Abou-Kassem (1975) model 


is only marginally better than the performance of the Dranchuk et al. (1974) model. 


One other mathematical representation of Standing and Katz's chart is a 
correlation developed by Gopal (1977). Although Gopal referenced more popular 
z-factor models, he presented few arguments to confirm the accuracy of his model. 
Because the accuracy of Gopal's model is unproven, it was not considered for this 


study. 


Since this study was concerned with the effect of acid gases and viscosity 
correlations on the development of Pseudo Pressure, only the original Dranchuk et al. 
(1974) EOS was used, even though the Starling EOS, the Hall and Yarborough EOS, 
and the Dranchuk and Abou-Kassem EOS are acceptable. 


The Dranchuk et al. (1974) EOS uses the Benedict, Webb and Rubin (BWR) 
EOS to develop the z-factor correlation. Although the BWR EOS can be quite 
awkward for general use in the field, the need for lengthy and repetitive calculations 
do not pose a problem for this study. In a situation in which it is not practical to use 
a multi-quadratic EOS, a simpler EOS which comprises a two constant cubic equation 


based upon van der Waal forces may be suitable. 


Models such as those presented by Peng and Robinson (1976), Schmidt and 
Wenzel (1980), and Soave (1972) improved the two constant cubic equation. While 
the EOS based upon van der Waal is more accurate than the BWR EOS in predicting 
PVT behaviour in the vapor/liquid equilibria, it is less accurate in predicting gas 
densities. Since this study dealt with single phase gas behaviour, the van der Waal 


EOS was not used. 
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1.2 GAS PROPERTIES AND CORRECTION FACTORS 


The most common and simplest mixing rule used to determine the critical 
properties for gas mixtures is Kay's Rule which considers the weighted average of 
the individual gas properties. Kay's Rule was used in this study with one exception. 
The exception to the use of Kay's rule is the use of Prausnitz and Gunn's rule of 
mixing in determining critical pressure. Dean and Stiel (1965) recommended the 
Prausnitz and Gunn method of evaluating critical pressure as it resulted in more 


accurate viscosity determination. 


From the critical properties, one can then determine the reduced properties 
which are in turn used to determine the z-factor, gas densities and compressibilities as 
long as the reduced values are within the model's limitation. However, if hydrogen- 
sulfide (H2S), carbon-dioxide (CO}), nitrogen (N2) or other non-hydrocarbon gases 
are present in the mixture, an additional correction will be required to match 
experimental z-factor, densities and compressibilities with the predicted results. 

Often correction for the presence of these gases is achieved by correcting either the 


z-factor itself or the properties which are used to predict the z-factor. 


Robinson et al.(1960), and Wichert and Aziz (1972) believed that in many 
situations only H2S and CO? gases are naturally present in sufficient quantity to 
affect gas behavior. The Robinson et al. method corrected the predicted 
compressibility factor by a ratio that was reflective of molar percent of CO? and HS. 
Later, Wichert and Aziz (1972) presented three mathematical correlations specific to 


the three popular z-factor models available at that time for these same gases. There 
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were, however, some limitations to both of these proposed methods of correcting for 
these acid gases. Firstly, the Robinson et al. (1960) charts are not suitable for 
computer applications, and if used they require a large datafile. Secondly, the 
Robinson et al. charts were also limited to a maximum combined CO? and H2S 
concentration of 60%. Finally, neither Robinson et al. or Wichert and Aziz presented 
a correlation to correct for the presence of N2 gas; in the case of Robinson et al., No 


gas was included with the methane. 


Despite its limitations, the Wichert and Aziz correlation (1972) was the best 
model to use in this study because it is applicable to computer programming and it 
provides a greater range of application for H2S and CO? gases than the Robinson et 


al. correlation. 


Although N? is naturally present in most reservoirs, the abnormally large 
concentrations which appeared in the sample may have been introduced by poor 
sampling techniques which allowed air to contaminate the sample. If contamination 
occurs, fluid analysis will indicate the presence of either O2 or oxidized 
hydrocarbons. In many cases, N2 represents a small proportion of the gas 
composition. However, to eliminate uncertainty, as well to eliminate the effect of this 
gas on the EOS, gas mixtures containing large amounts of nitrogen (greater than 


10%) were not used. 
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1.3 VISCOSITY OF NATURAL GASES 


Like compressibility, viscosity of natural gases has been studied by many 
investigators. The three models used in this study for an examination of viscosity 


correlations are, Carr et al. (1954), Dean and Stiel (1965), and Lee et al. (1966). 


Of these, the Carr et al. (1954) empirical correlation for viscosity appears to 
be the most popular in the Canadian petroleum industry. However, like the Standing 
and Katz z-factor chart, it is not readily adaptable to computer application. The 
Energy Resources Conservation Board (ERCB), however, has provided a program 
for the Carr et al. viscosity correlation in their third edition G-3 manual. The ERCB 
had also incorporated viscosity corrections for N2, CO2 and H2S gases, as originally 
presented in the Carr et al. (1954) paper. By comparison, no other viscosity models 


had a similar correlation for all three gases. 


Viscosity corrections for H2S and CO? in both the Lee et al. and the Dean and 
Stiel correlations were obtained through the use of the Wichert and Aziz(1972) 
correlation. The Wichert and Aziz correction factor was already incorporated in the 
correction to the gas properties, which had in turn corrected the EOS prediction. 
Both viscosity correlations require knowledge of gas density, which was calculated 


from values generated by the Dranchuk et al.(1974) EOS. 


Comparison of these three viscosity correlations was started by Chan and 


Dranchuk(1985). Although the results remain unpublished, the Chan and Dranchuk 
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study was regenerated in this study for the purpose of understanding the effects of 


different viscosity correlations in the treatment of acid gases. 


It was not the objective of the Chan and Dranchuk (1985) study nor was it 
the objective of this report to determine which of the three models is the most 
accurate. As well, the use of the Carr et al. (1954) viscosity correlation as the basis 
of comparison should not be construed as an indication that this correlation is the 
most accurate. The same caution is also expressed in the examination of the Pseudo 
Pressure determined with the use of the two viscosity correlations. In generating 
Pseudo Pressures, only the Carr et al. and the Dean and Stiel viscosity correlation 


were used due to high computer running cost. 


The Dean and Stiel (1965) viscosity correlation may have a greater range of 
application than Lee et al. (1966) because the Lee et al. study compared four gases in 
a limited number of temperature and pressure ranges. The Dean and Stiel correlation 
also has greater application than the Lee et al. viscosity correlation because Dean and 


Stiel considered both polar and non-polar gases. 


1.4 ANALYSIS OF THE DARCY-CONTINUITY (DIFFUSIVITY) EQUATION 


The three methods of gas transient analysis involve the Pressure, the Pressure 
Squared, and the Pseudo Pressure methods. As pointed out by Al-Hussainy and 
Ramey (1966), and later by Aziz et al. (1976), each method requires assumptions 


ohiied se Yo nninintiens ett ni eaten e qouke 207 strane seoen 
yolimeanss n) nnbiieanae wilesraern ia aku sa did boniiwessb surat 
AOTeLTTIOD eancealy (sit bri it ha ni iri od vine svc obese 
v9 ama Tuten dyin oF ah bau sa 


jo ager MieTg A svat yyiramnaiaienon vison ae!) sh? cannes? s 


nt eeey Wot Dewees bee lots 2a si apo 00!) Anieesd sodrnoliealigge 
halietamon lab? hag wor aT api ‘aurea bmmesningsrameon to secon bacon a 
ba cox sence noitelamado ytieooriy Lu ys, 990 sit suet notmOHgye eae ned ons 


vhs vekog- 00 bre raloq tied Gersbtanon fitz 


s 


= 
= 


about the reservoir properties. In all cases the reservoir is assumed to have the 


following properties: 


iv. 


iV. 


Vill. 


isotropic and homogeneous; 

isothermal; 

constant thickness throughout; 

horizontal; 

porosity, permeability and other formation properties are independent of 
pressure; 

gas with constant properties; 

laminar single phase fluid flows at a constant velocity; 

negligible effects from gravitational forces on the fluid; 

the formation fully penetrated by a vertical well; and 


the well open hole completed. 


Based on these simplifying assumptions, a second order, non-linear partial 


differential solution was derived as shown in Appendix I. However, even making all 


these assumptions, the equation cannot be solved in its present form. Using flow, 


boundary conditions and other simplifying assumptions, one can formulate the 


various Pressure, Pressure Squared and Pseudo Pressure solutions as indicated in 


Appendix I. Boundary conditions which are applied in the analysis of most PTA 


tests are: q = 


0 @ r=re and P =P; @t =0. In addition, the solution also is based on 


the assumption that the pressure transient acts as if it were in an infinite reservoir that 


has a wellbore of zero radius. 
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Homer PTA analysis, which accounts for the majority of analysis, is a two- 
rate test where the second rate is zero. The principle of super-position is used in this 


study to allow a multi-rate analysis. 


To solve gas PTA by the Pressure solution, it is assumed that the gas acts as a 
slightly compressible fluid (liquid). The same solution can be obtained if the product 


of P/uz is treated as a constant. 


In the Pressure Squared solution, it may be assumed that the gas is ideal, that 
viscosity is independent of pressure, and that the non-linear term of (dP/dr)? that is 
formed approaches zero rapidly; alternatively, the Pressure Squared solution can be 


solved by assuming that the product uz is constant. 


In tight or highly damaged reservoirs, (dP/dr)? may not approach zero rapidly 
and therefore should not be eliminated. This, however, means that the Pressure 
Squared solution cannot be solved because (dP/dr)? is a non-linear term. Dranchuk 
and Quon have proven that, in this situation, it is possible to solve the Pressure 
Squared equation without eliminating the (dP/dr)2 term, by defining new transform 


variables. 


Although these simplifying assumptions readily allow for the solution of the 
diffusivity equations, there is a question whether or not the assumption that P/j1z or 


1/1z are constant is correct, and, therefore, whether or not these simpler solutions 


should be applied. 
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1.5 PSEUDO GAS PRESSURE 


This section describes the development of the Pseudo-Pressure term. Since 
Al-Hussainy and Ramey (1966) first introduced their Pseudo-Pressure, there has 
been a debate as to whether or not the additional calculation, used to create greater 
accuracy, iS warranted in view of the greater simplicity of the Pressure or Pressure 
Squared method. Aziz et al. (1976) recommended which of the two methods would 


be the most suitable in specific situations. 


Aziz et al. (1976), indicated that by comparing the initial and final Pseudo 
Pressure with the appropriate Pressure or Pressure Squared curve, one can select the 
curve which best matches test response. If neither plot generates an appropriate 
curve for pressure range of interest, then, neither the Pressure nor Pressure Squared 


method will suffice. 


Dranchuk and Quon pointed out that the additional simplification of neglecting 
the non-linear term of (dP/dr)? in the Pressure Squared method was not necessary to 
solve the diffusivity equation. Large errors may be introduced in the Pressure 
Squared solution if the (dP/dr)2 term is inadvertently neglected. To avoid the 
problem, Dranchuk and Quon defined two transform functions to evaluate the 
(dP/dr)2 term. Depending on which variable is a strong function of pressure (either p 
or both p and 1) a specific transform function could be used. Dranchuk and Quon 
further indicated the error that could be expected, should this term be neglected. This 
error applies only to the differential equation and is not necessarily the same error 


which would appear in the final solution. 
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Zana and Thomas (1970) modified Pseudo Pressure by considering the 
effects of CO2 and H2S on real gas flows. The results were given as tables of 
Pseudo Pressure, with specified amounts of these gases. However, since Zana and 
Thomas used the Robinson et al. model to correct the gas property, their solutions 


were limited to a maximum concentration of CO? and H2S of 60%. 


Tiab (1984), tried to extend the effect of CO2 and H2S concentration by using 
the Wichert and Aziz correction for gas properties. Tiab, unfortunately, used the 
Gopal EOS to generate Z factors. The limitations of the Gopal EOS were discussed 


in Section 1.1. 


In 1987, Klins and Biterge used the following correlations to develop a semi- 
analytical solution to the Pseudo Gas equation: the Dranchuk et al. (1974) EOS, the 
Wichert and Aziz correction for CO? and H9S, and the Lee et al. correlation for 
viscosity. The resulting solution did not require numerical integration. For this 
study, similar attempts were made to develop a semi-analytical solution using the 
Dean and Stiel viscosity correlation, but these attempts were unsuccessful. In this 


study Pseudo Pressure was numerically integrated, using the trapezoidal rule. 


In all studies, including this study, the main objective is to increase accuracy 
in the determination and understanding of Pseudo Pressure in its application to PTA. 
Whether the Pseudo Pressure was determined by a semi-analytical solution or by 
numerical integration was not relevant to this study, provided that Pseudo Pressure 


was determined accurately. 


Although it is related to Pseudo Pressure, Pseudo-Time as defined by 


Agarwal, has not been considered in this study. Pseudo Time was first applied to aid 
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in the PTA analysis of hydraulic fractures. The application has since been extended 
to build up PTA of tight gas reservoirs and wells with changing storage coefficient. 
Since the wells which were analyzed were not tight and the time for which storage 


effects occur were limited, the use of pseudo time was not a concern in this study. 


As was mentioned earlier in this report, in PTA of gas reservoirs, three 
methods - Pressure, Pressure Squared, and Pseudo Pressure - are used in solving the 


Darcy-Continuity equation. 


Aziz, Mattar et al. (1976) compared the three approaches and determined the 
conditions under which the Pressure and Pressure Squared methods best approximate 
the Pseudo Pressure solution; they also determined the possible limitations of these 
two approaches. They also investigated the effects of using average properties on the 
resulting solution. However, neither Aziz et al. nor Al-Hussainy and Ramey, 
conducted any tests to determine whether or not the conditions which were assumed 
to be constant were actually constant. Nor did they discuss the effects either of acid 


gas or of the various viscosity correlations on the Pseudo Pressure determination. 


Aziz et al. (1976) used the Carr et al. (1954) model for viscosity in calculating 
the Pseudo Pressure. The Carr et al. model, however, is not easily adaptable to 
computer application as it requires storage of many data points. Also, some form of 
interpolation between points is required. The Energy Resource Conservation Board 
(ERCB) has developed and published a program that addresses these two problems. 
Firstly, the program enables one to use the Carr et al. graph in computer applications; 
secondly, the ERCB (1979) uses the Lagrange method for interpolation between data 
points. The Dean and Stiel (1965) model, however, is suitable for computer 


application because it is an algebraic expression. 
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CHAPTER II 


2.1 MODEL CONDITIONS 


In determining the various objectives of this study, computer programs for 
viscosity, z-factor, compressibility, density and Pseudo Pressure using Pseudo 


Reduced Pressures (pP;), were developed or obtained from available sources. 


Integration for Pseudo Pressure was accomplished by using the trapezoidal 
rule with a P; step interval of 0.05. Development of a semi-analytical solution for 
Pseudo Pressure, like the one presented by Klins and Biterge, was not possible due 


to the combination of empirical relations used in this study. 


The conditions were tested for a range of reduced pressures and temperatures 
for various gas mixtures and for actual pressure transient data available from the 
ERCB. The reduced temperature range varied from 1.05 to 3.00 with reduced 
pressure ranging from 0.20 to 20.00. Actual test data Pseudo Pressures were 


determined at the obtained pressure response and reservoir temperature. 


Thirty gas mixtures were considered for this study. They consisted of ten 
sweet dry gases, seven dry sour gases, four sweet wet gases, and nine sour wet 
mixtures. For the purpose of this study, wet gases were defined as those mixtures 
containing more than five mole percent of butane or heavier hydrocarbons. Although 


this may be less than an exact definition, these gases are most likely to produce 
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condensate. This definition provided a guideline permitting the classification 
of gases without knowledge of reservoir or surface pressures and temperatures. The 


gas composition of these mixtures is provided in Appendix II and at the beginning of 


each result output. 


Computer programs for the Carr et al. viscosity correlation and for the 
z-factor were provided in available programs in the Gas Well Testing Manual a 
edition and in a published paper by Dranchuk et al. (1974), respectively. In addition, 
the Carr et al. correction for No, H2S and CO3, was also accounted for in the 
viscosity correlation while the z-factor program was modified slightly to include 
compressibility determination. The other programs utilized were the Dean and Stiel 
viscosity correlation and the Wichert and Aziz correction for the acid gases H2S and 


CO>2. 


Programs used to determine P/j1z, [1z, cu, and Pseudo Pressure are 


provided in Appendix III, along with computer results. 


Some of the values generated by some of the programs were also plotted for 


comparison, using the University of Alberta Tell-a-Graf, Graphic Software Package. 


2.2 LIMITATIONS 


As stated earlier, results generated for all correlations used the same 


limitations as those established by the z-factor correlation. These limits were: 


1.05 < T; < 3.00; and 
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0.20 < P; < 20.0. 


For T; less than 1.05 and in particular in multi-phase regions, both the EOS 
and the Dean and Stiel viscosity correlation showed considerable deviation from 
actual behaviour. Problems can be minimized if the simpler cubic van der Waal EOS 
formula is used, as it predicts fluid behaviour more accurately in the multiphase 


region. 


Dean and Stiel recommend against using their viscosity correlation when the 
fluid has a reduced density greater than 2.0. This condition will occur at extremely 
low temperatures, or for gases within the two-phase region, or at extremely high 
pressure. In this study, high pressures, i.e., Py greater than 15, resulted in this limit 
being exceeded, but it was assumed that the relationship would be reasonably 
accurate for the purpose of comparison. Dean and Stiel had stated that "viscosity of a 


liquid at a higher reduced density depends on both density and temperature". 


For the Carr et al. viscosity correlation, similar limitations exist. These 
authors expressed concern about accuracy when T; approaches 1.0 and when P; is 
less than 1.0. In this region, the calculated ratio [/U1j is "expected to have the greatest 
uncertainty". Additionally, an upper limit of 860 degrees Rankin was set by the 
ERCB. In the analysis of the gas mixtures, this limit was exceeded, and again it was 
assumed that the correlation is reasonably accurate for the purpose of this study. 
Assumptions such as these are not uncommon as some well test analysis programs 


like those of Fekete make the same assumptions. 


The same reservation should be expressed in accepting available correction 


factors for Nz, CO, and HS as these correction factors are only approximate. Carr 
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et al. stated that such corrections only "indicate the direction and magnitude of error 


introduced and are hypothetical in nature...". 


In the Lee et al. viscosity correlation, no limits were expressed. However, 
since the Lee et al. study was based on the results of only four gas mixtures, through 
a temperature range of 100 to 340 degrees F and a pressure range of 100 to 8000 psi, 
the study was extremely limited. No attempts were made by Lee et al. to predict the 
viscosities of gases at extremely high pressure, or gas viscosity within the two-phase 
region. The authors indicated that the success of their viscosity prediction hinged on 
accurately predicting the gas density. It is assumed in the present study that gas 


density was accurately predicted. 
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CHAPTER III 


3.0 DISCUSSION 


Aziz et al. concluded that in certain conditions either the Pseudo Pressure 
solution or the Pressure solution can be used, as these solutions provide similar 
results: Aziz et al. also concluded that when certain other conditions prevail either the 
Pseudo Pressure solution or the Pressure Squared method can be used as these 
solutions provide similar results. Aziz et al. based their calculations on the 
assumption that the Pseudo Pressure solution is the most accurate. Mathematically 
and logically, the Pseudo Pressure solution is probably the most accurate method as it 
is the most consistent, that is, it requires the fewest assumptions. However, this 


does not guarantee that it will be the most accurate solution in all cases. 


It can be argued that some simplifying assumptions in the Pressure or 
Pressure Squared solutions may compensate for other parameters not accurately 
described in the solution. This argument is philosophical, with no one solution 
resolving the problem; it is best to study the actual field results to determine which 
theory is correct. In this study, as in the study conducted by Aziz el al., it was 


assumed that the Pseudo Pressure solution is the most accurate solution available. 


Examination of the products P/1z, 1/z and 1/cuL, was carried out on four gas 
mixtures, at four reduced temperatures. Each of these mixtures reflected extreme 


diversity in composition. The four mixtures were methane, a 53/47 percent split 
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methane/hydrogen sulphide mixture, mixture 19 (sometimes referred to as the 


wet mixture) with the following molar composition: 


GiGea67 718 n@d = 20:66 
ane ish iS ant 0x67 
Cage wel? :23 C6 - 0.41 
nGa> -_ 3.95 Ci aslo 
iC4 - 1.60 CBme- 0225 
eve O01 Noe) 0:60 
CQ? - 0.14 H2S - 0.00 


and mixture 23 (sometimes referred to as a wet sour mixture) with the 


following composition: 


Cl /'- 36.60 NCS wel LZ 
G2 =- 6.65 iC) 28 
(Ce aL C6 - 0.63 
nC4 - 4.44 Cree 053 
iC4 - 1.84 C8 - 0.00 
He - 0.00 lg eee CPA 
CO7 - 2.40 H5S. = 35.73 


While the gas products were plotted for the four gases using the Dean and 
Stiel viscosity correlation, viscosity and Pseudo Pressure were generated for 30 
different mixtures over similar pressure and temperature ranges, using two or more 
methods. The results of these calculations are presented in the Supplemntal volume 


of this study. 
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Gas mixture viscosities were determined by all three available correlations: 
Carr et al., Dean and Stiel, and Lee et al.. To examine the effects of the acid gases on 
viscosity, three conditions (treatments) were tested - using the available acid gas 
corrections, neglecting the acid gas calculated correction, and providing an analysis 


based on an acid free gas. 


Pseudo Pressures for the gas mixtures were generated using the Carr et al. 
and the Dean and Stiel viscosity correlations. The same three conditions that were 
tested in the determination of viscosity were also tested in the determination of 


Pseudo Pressure. 


3.1 ANALYSIS OF PRESSURE AND PRESSURE SQUARED ASSUMPTIONS 


The results of the P;/1z values calculated with the program in Appendix III 
were correlated on Cartesian co-ordinates as a function of reduced pressures (See Fig 
1). To investigate the trend of results, some of the same data was correlated on log- 
log co-ordinates (See Fig 2). Likewise, the inverse products of [tz were plotted 


against reduced pressure on Cartesian (Fig 3) and log-log (Fig 4) graphs. 


The results of this study indicate that the product, P;,/u1z, which is normally 
assumed to be constant (Appendix III), is not in fact constant. The results indicate 


that the Pressure solution method should not even be used in PTA of gas reservoirs. 


: 
a. 
>, a 
~ahotlalsrtad aldntuwe = ae 6 @ ilbihoon aisha ro ne sd 


ieee misy site aoutis st et leas ai ada si a — 


ton bron sitialiev ys oct gitar - bsieat sitter (ettionttnen) saation: vont ai 
a, 
éfeVieng fo aaihivotg brit ormeyrniog Sajal: any, boican ani guiosigas anol nexe) 


ney 998 bios caneue ite) 


oe Oe 3 


> v 


MAS Tia? Hi SAlsu Hoitisngy sigw aerikiiy 268 ala znuees ohveed 
si9W Ind) Sortnitbae) seul ae ST. Lanpdaleno. VNROOOIY isité bar enol od brs 
to momiienatsh oft nl Sota 'oes svew Vilgootw 6 nolieninneab oe Ai borat i 


sueeetl obese 


CYL 2G i AIO! Se JOR 9 Ovi 6 Si Jee 289 30 212 


iti ignites aq/. af riemers sri wy boralesiss novia aaa ool io Miya aT 


RPL942) esnucsare hess thea” ie bticr a ai “0 penile o9 near") w beaslorino ae . 


gal ne Hanalerrs anw cant owe ute arto onilsasn bron | agen a 
‘ 


brstealey ove? a2 Flo tetany; 3 Leh >it. Aefleadit an 208) eoonribn 8 got 


: — eT (6 919) gol-pal ban (fit) fiaieens2) rey cage nny vianiny 


w ren ial SAW ai henm  mabe a 0 atlases SAT 


pe ene patet at 3 Aa naam A) in 


a sae 


Even at high pressures P,/j1z is not constant, nor can it be approximated as a 
constant. Figure 1 indicate that when a combination of low temperatures and high 
pressures occurs, P;/|1z is the least sensitive to changes. By looking at the log-log 
graph in Figure 2, it can be seen that this product is not constant. Thus, even though 
P,/\1z is relatively insensitive to change at high pressure and low temperatures, it 


would be a mistake to approximate it as a constant. 


Figure 2 illustrates that the product P,/1z plots almost linearly against reduced 
pressure on both sides of the discontinuity, in a logarithmic relationship. As the 
temperature increases, the point of discontinuity also increases, so that ultimately the 


discontinuity disappears. 


When Aziz et al. applied the pressure solution, they determined the arithmetic 
average of pressures and then determined the gas properties using this arithmetic 
average. Use of this methodology requires caution because the assumed P,/1z 
product is not linear. For small pressure changes, a simple arithmetic average often 
suffices. However, for larger pressure ranges, the use of a logarithmic average 
(equation 3.1) is recommended. 


log fm) log m2) 


[Pr/mz] avg = 10 2 (350) 


The use of a logarithmic average is recommended, particularly at low 


temperatures, when P; is between 1.0 and 2.0. 


As in the Pressure solution, the Pressure Squared assumption with the 
product 1/1z being constant (Supplemental Volume), is equally incorrect. While the 
logarithmic plot of the 1/11z product (Fig4) appears to indicate that the product may be 


a constant, Cartesian plots (Fig 3) illustrate that this is not the case. The z-factor in 
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the product 1/\1z dominates this function at low pressures and low temperatures just 
as it did in the product P;/j1z. However, unlike the Pressure solution, the product 
1/\1z is less sensitive to temperature. As the pressure increases, the temperature 
differences between the solutions become less important. Therefore, where changes 
in temperatures would affect the Pressure solution, they would have less effect on the 


Pressure Squared solution. 


The usage of an average product would have greater application in this case 
for all conditions except one, namely, gases at low temperatures. Because the 
z-factor dominates the product 1/1z, only the logarithmic average should be applied at 


low temperatures. 
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3.2 ANALYSIS OF tc PRODUCT 


Compressibility is defined in equation 3.2. The domination of the z-factor in 
the calculation of compressibility (Equation 3.2), was first observed by Mattar et al. 
By using the Dranchuk et al. (1974) EOS, Mattar et al. were able to generate reduced 
compressibility versus a reduced pressure correlation. This plot differed from that 
previously published by Trube, in that it mathematically generated the slope dz/dP 
instead of relying on interpolation. Therefore, Mattar et a/. concluded that 
compressibility generated by the EOS should better predict the compressibility of 
fluids because compressibility generated by EOS is less subject to physical bias. 


dz 


dP 2) 


N |e 


el 
Ge 


However Figs 5 and 6 show that in the region where the function cyt appears 
to approach a constant, cy{t is not in fact a constant. Fig 6 illustrates the dominance 
of the z-factor in the product c;u. Mattar et al. reported this effect in their comparison 
with Trube's findings. In Fig 5, the Cartesian plot of this function does not show 


this problem due to the plot's scaling. 


In addition, as in the previous cases, a simple arithmetic average will often be 
sufficient to evaluate the gas property when small pressure ranges are investigated. 
However, if the pressure range is large, only a logarithmic average should be applied 
if accuracy is to be maintained. This recommendation does not apply for low 
temperature mixtures (T;<1.05) and reduced pressures between .8 and 2.0, due to the 


large discontinuity. In this range, actual PVT data or another EOS should be used in 
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order to accurately predict gas density and fluid viscosity. Then, a determination can 
be made as to which method, either an analytical method using another EOS, or 


actual data, will apply. 


Aziz et al., also indicate in their study when cy should be evaluated, i.e., at 
an average condition or initial condition. As discussed in the previous paragraph, the 
arithmetic average method will suffice when the range of pressure investigation is 


small, but the logarithmic average method would be more accurate. 


As demonstrated in Fig 7, if the pseudo gas equation is used, then, the 
product c;t appears to stabilize quickly to a constant. However, upon closer 
examination the ct product appears to reach a constant only at a low temperature 
and for a short pressure range. In addition, as in previous examinations, the Log- 
Log plot of the product cy, Fig 8, illustrates that the function still retains 
characteristics similar to those illustrated in Fig 6. Therefore, even though the 
function appears to stabilize quickly with the use of Pseudo Pressure, similar 


limitations will still apply as with the previous pressure vs. Cr plot. 


Although it is still common to find the use of the product cry in PTA, and the 
consequent difficulty of determining where to evaluate the product properly, it is 
recommended that Pseudo Time be used wherever possible. Agarwal, using the 
same logic as Al-Hussainy and Ramey had used in developing the concept of Pseudo 
Pressure, defined the integral Pseudo Time for gas PTA. Although not a true 
integral, the integral does account for changes in cyl product with respect to time at 
the corresponding pressure. This eliminates any difficulty in determining the point at 


which the product of cut should be evaluated. 
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However, as previously discussed, caution should still be exercised if the gas 
is examined within the two-phase region. Extension of the Mattar et al. 
compressibility study, Fig 6, illustrates that the observed inflection will become more 
predominant as the temperature decreases. When the gas enters the multi-phase 
region, the term dz/dP will be undefined, thereby dominating the compressibility 
function. Again, it is recommended in this situation that either actual PVT fluid 
analysis or compressibility determined by a different EOS, for example van der 


Waal's EOS, be used. 
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3.3 ANALYSIS OF H2S AND CO? ON THE GAS PRODUCT 


Several authors have attempted to correct the gas PVT behaviour for mixtures 
containing acid gases (H2S & COz) to those predicted by EOS. Two correction 
factors have been used in this study: the Wichert and Aziz acid gas correction, and 
Non-hydrocarbon gas correction within the Carr et al. viscosity correlation. 
Correction of the Dean and Stiel and the Lee et al. viscosity correlations, for the 
presence of acid gases, was achieved by correcting the gas density using the Wichert 


and Aziz correlation. 


CO and H2S gas treatment were handled in three ways - treatment "A", 
allowing for all correction factors; "B", eliminating these correction factors; and, "C", 
providing the analysis based on acid free basis. The type of method used referred to 
the viscosity correlation. Method 1 used the Carr et al. viscosity correlation, method 
2 used the Dean and Stiel viscosity correlation and method 3 used the Lee et al. 


viscosity correlation. 


In all cases, by examining the gas products P/uz, 1/\1z, and 1/cu, the function 
relation curve characteristic shape was retained throughout the complete pressure 
range at each specific temperature. Although the similarity of the curves is not 
apparent from an examination of the curves plotted in Fig. 1, 3, and 5, an 
examination of the Log-Log correlations (Fig. 2, 4, and 6), indicates that the curves 
do run parallel or nearly parallel to each other. The only exception to this rule was 
the plot of cyt product versus Pseudo Pressure (Fig 8). In this plot, the difference in 


gas composition results in a crossing of the curves at the low temperatures, rather 
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than a parallel plot. An examination of the effect of the gas correction may explain 


this phenomenon. 


If the gas composition is constant throughout the formation, the correction 
factors are constant regardless of changes to pressures or temperatures. For 
viscosity, as predicted by Carr et al. (1954), the gas correction factor is constant as 
shown in the ERCB program used in Appendix IV. However, the z-factor, 
compressibility and densities, used in the calculation of gas viscosities by Dean and 
Stiel and by Lee et al., is not corrected directly for H2S and CO>. Instead, the 
Wichert and Aziz correction factor, €3, which is a constant correction factor, is used 


to correct the critical properties used in the determination of c, 1, and p. 


When the Wichert and Aziz correction factor is applied, the critical properties 
of the mixture are shifted. This results in the shifts in the reduced pressures and 
reduced temperatures. In this study, with the exception of Fig 8, the reduced 
pressures and temperatures were held constant to show the effects of CO2 and H2S 


on viscosity correlation and Pseudo Pressure. 


For Pseudo Pressure, the results from method 1, using treatments A, B and 
C, and the results from method 2 using treatments A, B, and C, differ respectively by 
a constant. A similar relationship was observed in the changes in viscosity. This 
suggests the conclusion that, whatever correlation is used, when the data are 


processed in the same way, the results will differ only by a constant. 


However, as pointed out previously, and as indicated by Fig 8, actual field 
test data may not show this particular parallel response due to the presence of acid 


gases. When the reduced pressure and reduced temperature changes are due to the 
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shift in critical properties, the z-factor will be changed, thereby also affecting 
compressibility and density. Since the z-factor curves are not parallel, any shifts 
caused by the gas correction factor will not result in a parallel solution. At low 
temperatures and pressures, the effects of this shift will be the greatest, while at high 
temperatures and pressures such effects will be negligible. This can be seen in Fig 8 


> 


and is indicative of the result one expects to see in actual test results. 


3.4 EXAMINATION OF VISCOSITY CORRELATIONS 


This study is not a comprehensive study of the accuracy of available viscosity 
correlations, but rather an exploration of the effects of three treatments of acid gases 
on the viscosity results. The calculations of viscosity were made using viscosity 
correlations developed by Carr et al. (1954), Dean and Stiel (1965), and Lee et al. 
(1966). Correction for acid gases was carried out using the correction factors found 
in the Carr et al. viscosity correlation and by using the Wichert and Aziz z-factor 
correction which corrects gas density. Both Dean and Stiel and Lee et al. use density 


corrected by the Wichert and Aziz model. 


Although in this study the examination of the accuracy of the various models 
was not considered, a comparison of the methods was made to better understand 
differences in the results. The results of each calculation were compared with the 


results obtained using the Carr et al. method. The use of this method as the basis for 
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comparison should not be construed as an acceptance of the position that the Carr et 


al. model is more accurate than any other. 


All gas viscosity correlations react in a similar fashion to changes in 
pressures, temperatures, H2S and CO2. The overall trends are the same for all the 
models. As pressure increases, viscosity increases; as gas gravity increases, so does 
gas viscosity; as acidity of the gas increases, viscosity increases; as temperature 
increases, viscosity decreases. The rate of change to viscosity is also less at these 
higher temperatures. The exception to temperature changes occurred at low 
pressures, 1.e., less than 1.5: viscosity increases when reduced temperature 
increases. This is not surprising as this behaviour had also been shown in the 
correlations made by Carr et al.(1954) and by Lee et al. (1966). This information 
helps to explain the differences between the Pseudo Pressures arrived at when the 


various correlations are used. 


In this study, the different treatments of the acid gases resulted in parallel 
solutions when the Carr et al. viscosity correlations were used, and the same was true 
of the solutions obtained when the Dean and Stiel viscosity correlations were used. 
The Lee et al. method gave different results, in that the solutions were not parallel 


with each other. 


By treating H2S and CO? differently, i.e., by including all correction factors, 
by neglecting the correction factors, or by providing an analysis on an acid free basis, 
the viscosity results differed by a constant for all temperatures and pressures when 
the Dean and Stiel method was used. When the Carr et al. viscosity results for the 
three treatments were analyzed, the results differed by a constant for all pressure 


ranges as long as the temperature did not vary; there was a constant difference for 
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each temperature. Thus, while the constant was different at each temperature and for 
each gas mixture for Carr et al., the constant in the Dean and Stiel model remained the 
same for each gas mixture irrespective of the temperature. The Lee et al. results 
showed a different behaviour, viz., correction occurred in one direction only. The 


lowest difference, whether positive or negative, existed at the lowest pressure. 


An examination of the EOS explains the differences between the viscosities. 
Both the Dean and Stiel and the Lee et al. methods require knowledge of the gas 
density at a specific temperature and pressure. In this study density was determined 
by using the Dranchuk et al. (1974) EOS. Although knowledge of the gas gravity, 
reduced pressure and reduced temperature was necessary to calculate viscosity using 
the Carr et al. viscosity correlation, EOS was not necessary in the present study 
because the digitized tables of viscosity values provided by the ERCB were used. 
Interpolation between values in the table was achieved using the Lagrange 
interpolation method. However, as in the Trube presentation, much of the original 
data was smoothed to correct for any possible experimental errors and to present a 


graph appealing to the eye. 


By using the Dranchuk et al. EOS, both the Dean and Stiel and the Lee et al. 
methods showed a sinuousidal difference between their viscosity correlation and that 
of Carr et al.. In other words, the difference increases as reduced pressure increases, 


then suddenly decreases, only to increase again. 


In his unpublished work, Dranchuk notes that at low reduced temperatures, 
T,= 1.05, as reduced pressure increases, the slope dP/dV decreases to a minimum 


value, then starts to increase again. Since gas density is inversely proportional to the 


volumetric behaviour, a similar relationship between dP/dV and dP/dp holds true. 
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As the reduced temperature increases, the differences between the volumetric 
behaviours also decrease. That is, the effect of the multi-phase characteristic, the 
dew point or bubble point line, is reduced until a single changing dP/dV slope is 


observed. 


The increase in temperature has a similar effect when the differences between 
the viscosity results using the Carr et al. viscosity correlation, the Dean and Stiel 
viscosity correlation and the Lee et al. viscosity correlation are compared. As 
temperatures increase, the differences between the viscosity results from the three 
correlations decrease. However, the differences between the Carr et al. correlation 
and the Dean and Stiel correlation were generally smaller than the differences between 
the solutions reached using the Carr et al. and the Lee et al. methods. This was 
particularly the case for low temperature, high pressure gases containing large 
amounts of COy; in these cases, the differences between the Carr et al. and the Lee et 
al. methods were greater than 100 percent, as opposed to an 8 percent difference 
between the Carr et al. and the Dean et al. methods. The greatest deviation occurred 


near and above the Lee et al. pressure limit of 8000 psi. 
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3.5 COMPARISON OF THE RESULTS OF PSEUDO GAS PRESSURES 


In summary, comparison of the Pseudo Pressure results obtained by using 
the Carr et al. correlation with those obtained by using the Dean and Stiel correlation 
(Appendix V), results in the following conclusions. As temperature increases, the 
difference between Pseudo Pressures generally decreases, for sweet gases, and for 
gases that are wet or contain H2S. For gases that contain a high concentration of 


CO, the difference generally increases as temperature increases. 


In addition, the characteristics noted in the viscosity comparison are similar to 
the characteristics found in the Pseudo Pressure comparison. These characteristics 
are the constant difference due to the different treatment of acid gases and the 
sinuousidal fluctuation in the difference between the Pseudo Pressure using Carr et 
al. and the Dean and Stiel viscosity correlation. Considering the Dean and Stiel 


model only, the difference remained constant irrespective of the temperatures. 


In all cases, fluctuation in the difference occurred at low temperatures. 
However, as with viscosity, these fluctuations decreased as temperature increased. 
Even with mixtures containing high concentrations of CO2, such fluctuations were 
non-existent at the high temperatures. This indicates that, as with viscosity, the 


Pseudo Pressure determination is dominated by the EOS. 


Although the only variable in the determination of Pseudo Pressure was the 
viscosity model, the differences obtained differed from those obtained in the viscosity 


comparison. The difference between the Pseudo Pressures is not solely the result of 
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the inverse relationship of viscosity. Such inequality between the differences 
obtained in the viscosity correlation as opposed to the differences in Pseudo Pressure 


are the result of smoothing effects in the integration of the Pseudo Pressure. 


3.6 INTERPOLATION OF ACTUAL FIELD DATA 


As discussed earlier, the assumptions made within the Pressure and Pressure 
Squared methods were incorrect, and unless these solutions can be validated by 
actual field results, such methods should not be considered in PTA of gas reservoirs. 
Therefore, to compare results obtained by these methods to the Pseudo Pressure 
solution would be misleading and only confuse the issue. In addition, Aziz et al. 
(1976) have already shown the relationship between the Pressure, Pressure Squared 
and Pseudo Pressure solutions, with respect to the use of each solution. It is 


therefore not necessary to repeat the work in this study. 


Six actual field build up and deliverability tests were examined and the results 
are presented in Appendix VI. In all tests, Pseudo Pressures were generated using 
two viscosity methods as described in the Pseudo Pressure calculations. In the 
calculation of Pseudo Pressure, the pressure was fixed. This allowed pressure 
evaluation at different reduced pressures and temperatures. The reduced pressures 
and temperatures were different because three different acid gas treatments were 


used. Each gas treatment yielded slightly different critical properties. In summary, 
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the resulting reduced pressure and reduced temperature is different for each pressure 
and each gas treatment method. This allowed a full investigation of the effect of 


different gas treatments. 


Since the test data were supplied by the ERCB, it was not possible to 
scrutinize any information other than that actually provided by the testing company. 
Experience suggests that factors such as the incorrect reporting of fluid production, 
information as to well conditioning (particularly if this is a new completion), previous 
production history, and interference may influence test results. In addition, many 
tests required either a greater density of readings by rereading the copper strips or by 
using different and more accurate recorders such as electronic or quartz recorders. 
Since only a minimal number of readings were provided, it was impossible to fully 
analyse these tests. However, it was not the purpose of this study to reinterpret the 
build up tests, but rather to look at the effects of different gas treatments and viscosity 


correlations. Therefore, it was assumed that all tests were correctly reported. 


A typical pressure investigated during these tests was approximately 2 MPa 
(300psi), while the largest pressure was approximately 11 MPa (1550 psi). This 
large pressure change resulted in a reduced pressure range of approximately 5.5 - 
7.7. Although this range is limited, conclusions were based upon the results of the 


analysis. 


The Pseudo Pressures determined using the Dean and Stiel viscosity 
correlation (shown in the Supplemental Volume of this work in the Tabulated 
Results) resulted in consistently lower Pseudo Pressures than the Pseudo Pressures 
determined using the Carr et al. viscosity correlation. The differences between the 


two sets of results were relatively constant and were in all comparisons less than 0.25 
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of one percent. Even with the different treatments of the acid gases, the maximum 
deviation was less than 2.1 percent. However, unlike the results obtained using 
fixed reduced pressures, the results obtained using the fixed pressures showed that 
the differences increase as the pressure increases. Although these differences will 
also increase as the pressure range increases, the effects will not be noticeable on any 


PTA results. 


Since the PTA test results are dependent on the visual analysis of the Horner 
slope, such differences do not affect the accuracy of the analysis. If the analysis is 
limited to the region of transient pressure response, the reported range of differences 
is considerably reduced to 0.03 and 0.15 percent, respectively. However, during the 
early portion of the test, where PTA by Type Curve analysis was necessary, the 
difference may result in a slight decrease in the skin and storage coefficient. Because 
of the approximate nature of type curves, even with derivatives, and the number of 
curves available, the analysis will not differ even though the differences which occur 


in this region are greater. 


Although the Pseudo Pressures generated using different correlations differ, 
the solution for each test will be identical as long as the same Pseudo Pressure 


calculation method is used throughout the pressure range. 


Although the different treatments for acid gases yielded different Pseudo 
Pressures, resulting from the different reduced temperatures, the different treatments 
did not greatly affect the interpretation of the results. In the cases investigated, the 
differences between the methods were less than the differences due to the treatment of 


the acid gases. Where the previous section indicated that the treatment of acid gas 
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resulted in parallel Pseudo Pressures, the actual test data revealed that the differences 


increased slightly as pressure increased. 
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CHAPTER IV 


4.0 CONCLUSION 


The purposes of this study were to examine the assumption that P/1z is 
constant in the pressure equation and the assumption that 1/uz is constant in the 
Pressure Squared equation. The results of the study showed that these assumptions 
were incorrect. If either the Pressure or Pressure Squared equation must be used, 
then, a logarithmic average pressure should be used to accurately evaluate the gas 
parameters, i.e., viscosity, compressibility, and the compressibility factor. A simple 
arithmetic average pressure should not be used unless the pressure range of 


investigation is relatively small. 


Although the treatment of acid gases appeared to have little effect on the 
outcome of the solution, this does not mean that correcting gas properties for acid 
gases should never be considered. In areas other than PTA, correction for these 
gases to actual PVT behaviour is still necessary. In this study, however, it was not a 
concern because PTA is concerned with the relationship between pressure differences 


and not with the pressure itself. 


Equally, this study showed that the choice of EOS methods by which Pseudo 
Pressure is determined has no effect on the PTA solution for the same reasons. 


However, as indicated above, not all EOS used in PTA will produce parallel 
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solutions. Limits of the particular EOS may restrict the use of some 
correlations, for example in the multi-phase region. This study did not address this 


problem. 


The results of the three viscosity correlations differed from each other and the 
results of the two Pseudo Pressure calculations also differed from each other: these 
calculations were made using a constant reduced temperature. The magnitudes of the 
differences were not compared as they varied with the gas composition (refer to 
Appendixes IV and V). Often, the differences between the methods were constant at 
high temperatures, regardless of the different treatments for acid gases. This resulted 
in parallel solutions. At low temperatures, there was fluctuation in the differences 
between the viscosity correlations because both the Dean and Stiel method and the 
Lee et al. method depend on the EOS in the calculation of viscosity. The differences 


in the Pseudo Pressure results showed the same type of fluctuation. 


When actual field results were analysed, reduced temperatures were different 
because of the different critical temperatures due to the type of acid gas treatment. 
The results, although slightly different, also indicated a relative consistency between 
methods over the pressure range of investigation. These results translated to parallel 
solutions, and therefore the same answer, even though the difference increased 
slightly as pressure increases. The three methods of treating the acid gases resulted 
in greater differences than the differences between the correlations themselves, but 
these difference were still relatively constant. The difference in the constant is 


unlikely to noticeably affect any PTA results. 
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PTA interpretation of actual field tests is limited by two factors, namely, the 
quality of the data and the skill of the interpreter. Often, data quality is dependent 
upon the accuracy of the recording device: accuracy may range from 0.5% FSD for 
mechanical (Amerada) recorders, to 0.005-0.05% FSD for the electronic recorders. 
In addition, data analysis still requires human interpolation, which in itself may 
account for greater differences than the differences between the methods or in the 
treatment in acid gases. Although the advent of automatic curve fitting software and 
PTA programs may ultimately eliminate the human factor, perfection of this software 
is still some time away, even though several companies are making remarkable 


progress. 
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CHAPTER V 


5.0 RECOMMENDATIONS 


It is recommended that data from several fields be analysed, each reflecting a 
completely different fluid composition. In selecting the particular tests, consideration 
should be given to selecting tests with three distinct flow regimes - storage, transition 
and transient. Further analysis of these tests will likely confirm the observations in 


this study. 


It is recommended that the Carr et al. (1954) correction factor be investigated 
to confirm its accuracy, and that its results be compared with the results obtained 


using Wichert and Aziz correction factor. 


It is also recommended that a comparison of all viscosity correlations be made 
to verify the accuracy of existing correlations. Although the use of the Dean and Stiel 
and the Carr et al. viscosity correlations produced little difference in the results of the 
PTA analysis, as this study has shown there can be extreme differences between the 
results of these correlations when there are large concentrations of acid gases. These 
differences can affect metering factors, pipe designs, and fluid flow characteristics, 


which depend on accurate values. 
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APPENDIX I 


DERIVATION OF THE PRESSURE AND PRESSURE SQUARED EQUATION, 


Fluid flow in a reservoir has been described for many years by the Darcy- 
Continuity equation. In developing this equation the following conditions had been 
looked at. These are: 

i, a constitutive equation which describe and relates the shear stress imposed 
on the fluid to the shear rate; as with the equation which describe the 
viscosity of a Newtonian fluid at any temperature or pressure: 

ii. the force momentum balance of the system, as described by Newton's 
second law of motion: 

ili. the conservation of mass: 

iv. the equation of state which relates fluid density to temperature and 
pressure: 


v. and the conservation of energy. 


In addition to these assumptions, it was also assumed that the reservoir or 
fluid would be: 

i. Single phase fluid with constant property; 

ii. formation which is homogeneous, isotropic and incompressible; 

iii. isothermal and laminar fluid flow; 


iv. reservoir of uniform thickness; 
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Vv. premeability and porosity independant of pressure; 
vi. no cross flow between beds or outside top and bottom bounderies; 


Vii. ‘gravitational effect which are minimal. 


Considering the radial case only, and the previous assumption, one can derive 


the Darcy Continuity equation as expresed below; 


v: Ra =-( 52) (1) 


where, 


E=g-5 VP (2) 


By subsituting equation 2 into 1 and expressing it in terms of polar cordinates 


the following equation can be derived; 


rue Tpke east 
rile @)" Ce) ° 


However, for gas, density can be described as; 


Therefore; 
1 0 /Pk @® dP 
se del gg AS, ES oa Aad (5) 
Lor par ar | ; raz) 


Klinkenberg however proved that permeability was not constant but varies 
inversly with pressure. He had expressed permeability of gas flow through a porous 


media as; 
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Aronofsky and Ferris later indicated that the equation 5 depended more on 


changes in gas properties than changes in premeability. Therefore Klinkenberg 


effects were neglected in this study. 


By simplifying further, equation 5 becomes; 
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CASE 1 PRESSURE METHOD 


In this case, gas which is the only fluid present, is treated as a liquid. Two 
assumption can be made to simplify the equation. First the product of P/\z is 
assumed to be constant, or second the pressure gradient, (dP/dr)2, is assumed to be 
small. If these were the case equation 6 may be reduced to the same result. The only 


simplified presented is for case one as follows; 
ge <0P 0 /P 
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Since the gas is treated as a liquid or a slightly compressible fluid, Z Factor is 


assumed to be constant, with compressibility being expressed as; 


1 
a, (8) 
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Therefore equation 7 becomes; 
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By expanding equation 9 to get 


027P 1 0P uc OP 
— += — = 10 
o2 OT Or)hCUkK Ot (10) 


CASE 2 PRESSURE SQUARED METHOD 


In this method, gas is treated ideally. With this method either the product of 
1/uz is assumed to be constant, or that 1 varies independantly with pressure since 


z= 1 (ideal). Knowing this equation 6 can be simplified to: 
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Letting 

d(P?) = 2PdP (12) 
Equation 11 becomes; 

1 9 (r H(P2)) _ op o(P?) (13) 
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Again equation 8 still applies so that equation 13 becomes 
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Expanding equation 14: 
d2(P2) 1 a(P2 oP o(P2 
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Assuming the pressure gradient (dP/dr)” is small and therefore neglectable 


then equation 15 becomes; 
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CASE 3 PSEUDO-PRESSURE 
By using Al-Hussainy et al. definition of pseudo-pressure 
Pap 
yO) | = ap (17) 
ak 
one can redefine equation 6 as 
1 0 (r OF) _ o OP/z) (18) 
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and compressibility expressed by 
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So that substituting this into equation 20; 


Expanding equation 23 
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APPENDIX III 


THIS PROGRAM IS USED TO GENERATE VALUES FOR 
-REDUCED COMPRESSIBILITY 
-COMPRESSIBILITY FACTOR 
- VISCOSITY 
-PRODUCT OF COMPRESSIBLITY AND VISCOSITY 
-PSEUDO GAS PRESSURE 


RESULTS ARE ATTACHED 


METHOD OF NUMERICAL INTERGRATION IS BY SIMPSON"S 
1/3 RULE FOR EVEN STEP SIZES AND A COMBINATION OF 
SIMPSON'S 1/3 & 3/8 RULE FOR ODD STEP SIZES. 


VISCOSITY HAS BEEN CALCULATED USING D.DEAN & L.STIEL 
MODEL. THE PAPER WAS PRESENTED IN AIChE JOURNAL, 
MAY 1965. 


Z FACTOR, REDUCED COMPRESSIBILITY AND DENSITY HAS 
BEEN CALCULATED USING DRANCHUK, PURVIS, & ROBINSON 
MODEL WHICH WAS PRESENTED IN THE TRAN OF AIME VOL. 146, 
PP140-149 


SOUR GAS CORRECTION BY WICHERT & AZIZ CORRELATION 


DOUBLE PRECISION NTC,NPC,SM,CG,XX 
DOUBLE PRECISION XPROD,PROD,VG,VISC 
DIMENSION P(20) 


READ (8,10) G,PC,TC,T 


10 FORMAT (F7.5,3(2X,F8.4)) 


SM = G*28.75 
READ (8,100) (P(1J), IJ=1,19) 


100 FORMAT (F5.1,1X,9(F6.1,1X),/,9(F6.1,1X)) 


1 


TR =T/TC 
L=2 
H=0.05 
WRITE (5,1) TR 
WRITE (5,50) 
WRITE (5,2) 
WRITE (5,3) 
FORMAT ('1'//,2X,'The Value of Tr is = ',F5.2) 


50 FORMAT (//) 


2 


FORMAT (/,10X,'Pr',8X,'Pr/CrU',7X,'U',6X, N/GrU GX, 


63 


a V2 TMT OR Ch Vig y vate 7s 
MC) 4 j 101% PAWIE ‘ sax IDA Sis 4 
ZANI2 WTS Oe 


Te. PAB 
AAs JAN At 


BAHSCTIRVISC CsA Y TLL eae 
AOZALAOR D 21VHIY AAA Aer 
abl. JOV SMIA 90 MAST UAT y: 


CVO 


OOM 


64 


TOL Oren) VISCO 5), Cr',7X, Pr/ZU',6X, 
TZU oe. Dr) 
FORMAT (4X,106(‘_'),/) 


NTC=TC*5.0/9.0 

NPC=PC/14.7 

ZI=NT C**(1.0/6.0)/(SM**0.5*NPC**(2.0/3.0)) 
IF (TR-1.5) 5,6,6 

VR=34.0E-05*TR**(8.0/9.0) 

GO TO7 
VR=166.8E-05*((0.1338*TR-0.0932)**(5.0/9.0)) 
PR=0.2 


DO 4 I=1,19 
PX=P(I) 

PR=PX/PC 
N=((PR-0.195)/H)+2 
IF (N.LE.2) GO TO 33 
J=N-(N/2)*2 


IF (PR.GT.30.0) GO TO 4 
H=(PR-0.2)/N 

SUM=0.0 

XINTG=0.0 

PR1=0.2 


IF (J.EQ.0) GOTO 8 
SIMPSON'S 3/8 RULE 


CALL ZKATZ (PR1,TR,DR,ITER,Z,CR) 
CALL VISC (ZI,VR,DR,V) 
SUM=PR1/(V*Z) 

DO 11 =1,3 

PR1=PR1+H 

CALL ZKATZ (PR1,TR,DR,ITER,Z,CR) 
B=1.0 

IF (ITER.EQ.0) GO TO 24 

CALL VISC (ZI, VR,DR, V) 

IF (11.EQ.3) GO TO 9 

A=3.0 

GO TO 10 


9 A=1.0 
10 SUM =SUM + A*(PR1/(V*Z)) 


11 


CONTINUE 
XINTG=3.0*H*SUM/8.0 
L=4 
N=N-1 
IF (N.LT.L) GO TO 34 
SUM=0.0 
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C SIMPSON'S 1/3 RULE 
o 


8 CALL ZKATZ (PR1,TR,DR,ITER,Z,CR) 
B=2.0 
IF (ITER.EQ.0) GO TO 24 
CALL VISC (ZI, VR,DR,V) 
SUM = PR1/(Z*V) 
DO 12 KK=L,N 
PR1=PR1+H 
CALL ZKATZ (PR1,TR,DR,ITER,Z,CR) 
CALL VISC (ZI, VR,DR,V) 
M=KK-(KK/2)*2 
IF(M.EQ.0) GO TO 13 
A=2.0 
GO TO 14 
13. A=4.0 
14. SUM=SUM+A*(PR1/(Z*V)) 
12 CONTINUE 


PR1=PR1+H 
CALL ZKATZ (PR1,TR,DR,ITER,Z,CR) 
B=4.0 

IF (ITER.EQ.0) GO TO 24 

CALL VISC (ZI, VR,DR,V) 
SUM=SUM+PR1/(V*Z) 


XINTG=XINTG+H*SUM/3.0 


GO TO 34 
33. XINTG=0.0 


34. CALL ZKATZ (PR,TR,DR,ITER,Z,CR) 
B=5.0 
IF (ITER.EQ.0) GO TO 24 
CALL VISC (ZI,VR,DR,V) 
U=XINTG 
CRU = V*CR 
PRCRU=PR/CRU 
XCRU=1.0/CRU 
CG=CR/PC 
PRZU=PR/(Z*V) 
XZU=1/(Z*V) 
WRITE (5,15) PR,PRCRU,U,XCRU,CRU,V,CR,PRZU,XZU,DR,PX, 


+ Z,CG 
15 FORMAT (' ' 2X,F10.2,2X,6F10.4,2X,4F10.4,2F10.7) 


GO TO 4 
24 WRITE (7,25 | 
25 FORMAT (T ‘ #4 TR OR PR EXCEEDS MODEL LIMITS. ***") 
WRITE (7,26) PR,TR,B mee 
26 SORMAT (2X,'PR=',F7.4,6X,'TR=',F6.4,6X,'B= ',F4.1) 
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GO TO 40 


4 CONTINUE 


2 


CONTINUE 


40 STOP 


os 
C 


ele 


END 
SUBROUTINE VISC (ZI, VR,DR,V) 


VG=10.8E-05*(EXP(1.439*DR)-EXP(-1.11*DR**1.858)) 
V=(VR+VG)/ZI 

RETURN 

END 


SUBROUTINE ZKATZ(PR,TR,DR1,ITER,Z,CR) 


DIMENSION A(8) 

DATA A/0.31506237,-1.0467099,-0.57832729,0.53530771, 
= -0.6123032,-0.10488813,0.68157001,0.68446549/ 
ITER=0 
DR=1.0 

IF(TR-1.05)10,1,1 


1 IF(TR-3.0)2,2,10 
2 IF(PR-30.0)3,3,10 
3 DO 8 ITER=1,10 


DR2=DR**2 

T 1=(A(1)*TR+A(2)+A(3)/TR**2)*DR 
T2=(A(4)*TR+A(5))*DR2 
T3=A(5)*A(6)*DR**5 

T4=A(7)*DR2/TR**2 

T5=A(8)*DR2 

T6=EXP(-T5) 

P=(TR+T 1+T2+T3)*DR+T4*DR*(1.0+T5)*T6 
DP=TR+2.0*T 1+3.0*T2+6.0*T3+T4*T6*(3.0+3.0*T5 
+ -2.0*T5**2) 

DR1=DR-(P-0.270*PR)/DP 

IF(DR1)4,4,5 


4 DR1=0.5*DR 
5 IF(DR1-2.2)7,7,6 

6 DR1=DR+0.9*(2.2-DR) 

7 IF(ABS(DR-DR1)-0.1E-04)9,8,8 
8 DR=DR1 

9 Z=0.270*PR/(DR1*TR) 


DZ = (T1+2.0*T2+5.0*T3+2.0*T4*(1.0+T5 
+  -T5**2.0)*T6)/(TR*DR1) 
CR = 1.0/PR -0.270*DZ/(Z**2*TR*(1.0 + DR1*DZ/Z)) 


10 RETURN 


END 
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APPENDIX IV 


THIS PROGRAM IS USE TO DETERMINE VISCOSITY 
BY DEAN AND STIEL MODEL, CARR, KOBAYASHI 
AND BURROWS, AND LEE GONZALEZ AND AND EAKIN; 


REDUCED DENSITY IS SUPPLIED BY 
DRANCHUK, PURVIS, AND ROBERTSON'S 
MODEL OF THE BWR SOLUTION. 


THIS PROGRAM ELIMINATED ALL EFFECTS OF CO2 ,H2S AND N2, 
FROM THE GAS COMPOSITION. 
METHOD A INCLUDES ALL CORRECTION FACTORS, 
METHOD B EXCLUDES THESE CORRECTION FACTORS, 

AND METHOD C PROVIDES ANALYSIS ON ACID GAS FREE BASES. 


CARR, KOBAYASHI AND BURROWS CORRELATION AND LAGRANGE 
INTERPOLATION WERE PROVIDED BY THE ERCB IN G-3 
APPENDIX 'D', 3 RD EDITIONS 


RESULTS OF THIS PROGRAM IS ATTACHED 
NOMENCLATURE USED WITHIN THIS PROGRAM 


A - MOLAR FRACTION OF H2S + CO2 

B - MOLAR FRACTION OF H2S 

DE - GAS DENSITIES 

DIFF- DIFFERENCE BETWEEN THE VARIOUS METHODS 

E3 - WICHERT & AZIZ CORRECTION FOR ACID GASES 

G - GAS GRAVITY (ENTER 0.0 IN DATA READ STATEMENT 
TO CALCULATED G FROM COMPOSITION) 

NN - NUMBER OF MIXTURES 

NP - CRITICAL PRESSURE ATM 

NT - CRITICAL TEMPERATURE K 

P -EVALUATION PRESSURE PSIA 

PPC - CRITICAL PRERRURE PSIA 

PR - REDUCED PRESSURE 

R_ - UNIVERSAL GAS CONSTANT 
10.732 PSIA-CU-FT/DEG-R-LBMOL 

SM - PSEUDO MOLECULAR WEIGHT LB/LBMOL 

TC - CRITICAL TEMPERATURE DEGR 

TCO - OLD CRITICAL TEMPERATURE FOR USE IN WICHERT & 
AZIZ CORRELATION 

TR - REDUCED TEMPERATURE 

VC - CRITICAL VOLUME CU-FT/LBMOL 

V1 -GAS VISCOSITY AS CALCULATED BY DEAN & STIEL 


METHOD; cp 
V2 -GAS VISCOSITY AS CALCULATED BY CARR, KOBAYASHI 


67 


Lvs @f, SOD 40 2 Tose 


yeh MOr Oe 
Pa Bier sn 


292A8 Sat SAA ag AEA 


fA RAL CAA MOET AIIAROR i 
_-OViSO89 ST VA oa 


CHROATTA nat t 
MARIONT ill 


49 wa 
a 


‘ Th ie ao ar 


2721-4 ee 


: artored >: : 


PlSpeleselelolotclelolelolelole(elele lo lel orelele lele lee lee eloiGlere @Krel@n-(ele 


68 


& BURROWS; cp 
V3 - GAS VISCOSITY AS CALCULATED BY LEE, GONZALES, 
& EAKIN; cp 
Rone mae RATIO FROM CARR, KOBAYASHI & BURROWS 
VG - VISCOSITY FACTOR LESS THE VISCOSITY EFFECTS AT 
MODERATE PRESSURE (0.2-5 ATM) 
VR - RESIDUAL VISCOSITY 
WM - MOLECULAR WEIGHT LB/LBMOL 
X(1)-X(14)-MOLAR % OF GAS COMPOSITION IN THE ORDER OF 
C1,C2,C3.1C4,NC4,IC5,NC5,C6,C7,C8,CO2,H2S,N2,HE 
XK, - CONSTANTS FOR LEE, GONZALEZ AND EAKIN 
XX, WISCOSITY CORRELATION 
DS & 
XT - SANDFACE TEMPERATURE DEG R 
ZC - CRITICAL COMPRESSIBILITY FACTOR 
ZI - VISCOSITY PARAMETER 
ZP - PRESSURE, PSIA 


PROPERTIES LISTED BELOW HAVE BEEN OBTAINED FROM 
GAS PROCESSOR SUPPLIERS ASSOCIATION, ENGINEERING 
DATA BOOK, NINTH EDITION 


CONVERSION USED TO OBTAINED UNITS USED IN THE PROGRAM 
ARE: 

TC = TC(OF) + 460.0 

VC = VC(FT3/LB)*WM(LB/LBMOL) 

ZC=/PG*EVC/R*TC 


GAS COMPOSITION IF KNOWN, SHOULD BE EXPRESSED IN 
PERCENTAGES AND ENTERED IN THE FOLLOWING ORDER 


TAC DAC3, nC4, iC4, nC5, 1C5, C6, C7 
C8+, CO2, H2S, N2, He 


CRITICAL PROPERTIES HAVE BEEN DETERMINED BY KAY'S RULE 
EXCEPT FOR CRITICAL PRESSURE WHICH HAS BEEN 
DETERMINED BY PRAUSNITZ AND GUNN MIXING RULE. 


DIMENSION ZPC(5),ZTC(5),ZERR(5),ZF(5),ZT (5), TT (5) 
DIMENSION XP(5),X(14),TP(5),ZP(6),F(6) 

DOUBLE PRECISION ZX1,ZX2,ZX3,VC,SM,SVC,TT,XP,V1,V2,Z 
DOUBLE PRECISION G1,G2 


READ (5,5) SSA,NN 
FORMAT (F6.3,1X,12) 
TR=1.05 

R=10.732 
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ATURE AR YR Ce . 
ALIA 3 


WRITE (6,10) TR 
WRITE (7,11) TR 

10 FORMAT (2X,'THE FOLLOWING CALCULATIONS ARE FOR A ' 
+'REDUCED TEMPERATURE OF ',F4.2,//) 

11 FORMAT (2X,'THE FOLLOWING CALCULATIONS ARE FOR A ' 
+'REDUCED TEMPERATURE OF ',F4.2,//) 

CG 


> 


> 


DO 240 I=1,NN 

20 READ (5,30) (X(1J), U=1,14) 

30 FORMAT (14(F5.2,1X)) 

fe 
CALL MIXT(X(1),X(2),X(3),X(4),X(5),X(6),X(7),X(8) 
+,X(9),X(10),X(11),X(12),X(13),X(14),ZPC(1),ZT C(1), 
+H2S1,CO21,XN21,ZPC(2),ZTC(2),H282,CO22,XN22, 
+ZPC(3),ZTC(3),XG1,XG2,SM1,SM2) 


C 
C CHECKING UPPER AND LOWER BOUNDARY AS FOLLOWS 
S 


DO 35 IA=1,3 
ZERR(IA)=0.0 
TR=1.05 
CALL BOUND (TR,ZTC(IA),XG1,XG2,ZERR(IA),ZT (IA)) 
35 CONTINUE 
fe 


PR=0.25 


C CONSTANT FACTOR DETERMINATION FOR DEAN AND STIEL 
C NATUAL GAS VISCOSITY CORRELATION 
G 
IF(TR-1.5) 80,90,90 
80 VR=34.0E-05*TR**(8.0/9.0) 
GO TO 100 
90 VR=166.8E-05*(0.1338*TR-0.0932)** (5.0/9.0) 
100 DO 105 IC=1,3 
TT (IC)=ZTC(IC)*5.0/9.0 
XP(IC)=ZPC(IC)/14.7 
105 CONTINUE 
Z1=TT(1)**(1.0/6.0)/(SM1**0.5*XP(1)**(2.0/3.0)) 
Z2=TT (2)**(1.0/6.0)/(SM1**0.5*XP(2)**(2.0/3.0)) 
Z3=TT (3)**(1.0/6.0)/(SM2**0.5*XP(3)**(2.0/3.0)) 


Ova pl 5) Oy wh Ve os 2a) 2 noni 
C CONSTANT FACTORS DETERMINATION FOR LEE GONZALEZ AND 


EAKIN 
C NATURAL GAS VISCOSITY COORELATION AS PRESENTED IN SPE AUG 
C 1966. 
G 
XT1 = TR*ZTC(1) 
XK1 = ((9.4 + 0.02*SM1)*XT1**1.5)/(209.0 + 
+ 19.0*SM1+XT1) 
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XX1 = 3.5 + 986.0/XT1 + 0.01*SM1 
XY1 = 2.4 - 0.2*XX1 


XT2 = TR*ZTC(2) 

XK2 = ((9.4 + 0.02*SM1)*XT2**1.5)/(209.0 + 
+ 19.0*SM1+XT2) 

XX2 = 3.5 + 986.0/XT2 + 0.01*SM1 

XY2 = 2.4 - 0.2*XX2 


XT3 = TR*ZTC(3) 

XK3 = ((9.4 + 0.02*SM2)*XT3**1.5)/(209.0 + 
+ 19.0*SM2+XT3) 

XX3 = 3.5 + 986.0/XT3 + 0.01*SM2 

XY3 = 2.4 - 0.2*XX3 


DO 205 IB=1,50 
Te(PReet OOF) GOTO 160 
I (PR..Gi620,0):GO.TO 220 


XP1 = PR*ZPC(1) 

XP2 = PR*ZPC(2) 

XP3 = PR*ZPC(3) 

XT1 = TR*ZTC(1) 

XT2 = TR*ZTC(2) 

XT3 = TR*ZTC(3) 

CALL ZKATZ (PR,TR,ZZ,DR,ITER) 

DE1 = XP1*SM1/(R*62.4*XT 1*ZZ) 

DE2 = XP2*SM1/(R*62.4*XT2*ZZ) 

DE3 = XP3*SM1/(R*62.4*XT3*ZZ) 

CALL VIS3 (XK1,XX1,XY1,DE1,V3A) 

CALL VIS3 (XK2,XX2,XY2,DE2,V3B) 

CALL VIS3 (XK3,XX3,XY3,DE3,V3C) 

CALL VIS2 (PR,TR,XG1,ZT(1),H2S1,CO21,XN21,V1A) 
CALL VIS2 (PR,TR,XG1,ZT(2),H2S2,CO22,XN22,V1B) 
CALL VIS2 (PR,TR,XG2,ZT(3),H2S2,CO22,XN22,V1C) 
CALL VIS1 (Z1,VR,DR,V2A) 

CALL VIS1 (Z2,VR,DR,V2B) 

CALL VIS1 (Z3,VR,DR,V2C) 


DIFF1 = (V1A-V1B)*100.0/V1A 
DIFF1A= (V1A-V1C)*100.0/V1A 
DIFF2 = (V1A-V2A)*100.0/V1A 

DIFF2A= (V2A-V2B)*100.0/V2A 
DIFF2B= (V2A-V2C)*100.0/V2A 
DIFF3 = (V1A-V3A)*100.0/V1A 

DIFF3A= (V2A-V3A)*100.0/V2A 
DIFF3B= (V3A-V3B)*100.0/V3A 
DIFF3C= (V1B-V3B)*100.0/V1B 
DIFF3D= (V2B-V3B)*100.0/V2B 
DIFF3E= (V3A-V3C)*100.0/V3A 
DIFF3F= (V3B-V3C)*100.0/V3B 
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DIFF3G= (V1C-V3C)*100.0/V1C 
DIFF3H= (V2C-V3C)*100.0/V2C 
Cc 
DO 140 IC=1,3 
ZP(IC) = PR*ZPC(IC) 
140 CONTINUE 
G 
WRITE(6,155) LPR,ZP(1),ZP(2),ZP(3),V1A,V1B,V1C, 
+V2A,V2B,V2C,V3A,V3B,V3C 
CALL ZKATZ (PR,TR,ZZ,DR, ITER) 
C 
WRITE(7,156) I,PR,DIFF1,DIFF1A,DIFF2,DIFF2A,DIFF2B 
+, DIFF3, DIFF3A ,DIFF3B,DIFF3C,DIFF3D,DIFF3E,DIFF3F, 
+DIFF3G,DIFF3H 
155 FORMAT(I2,2X,F5.2,2X,17(F10.4,2X)) 
156 FORMAT(I2,2X,F5.2,2X,17(F8.4,2X)) 
G 


IF (PR-0.29) 160,165,165 
160 PR=PR+0.05 

GO TO 205 
165 IF (PR-1.49) 170,175,175 
170 PR=PR+0.10 

H=0.02 

GO TO 205 
175 IF (PR-2.99) 180,185,185 
180 PR=PR+0.50 

H=0.1 

GO TO 205 
185 IF (PR-9.99) 190,195,195 
190 PR=PR+1.0 

H=0.20 

GO TO 205 
195 IF(PR-19.999) 200,240,240 
200 PR=PR+5.0 

H=0.5 
205 CONTINUE 
220 WRITE (6,212) 

WRITE (7,212) 
212 FORMAT (''/) 
240 CONTINUE 


G 
WRITE (6,211) ZTC(1),ZT C(2),ZTC(3),ZPC(1),ZPC(2), 
+ZPC(3) 
WRITE (7,210) XT1 . 
210 FORMAT (//,, TEMPERATURE IS OUTSIDE RANGE OF CKB’, 
+' MODEL WHERE TEMP>860:',/, IT IS ASSUMED THAT 
+ VISCOSITY',' CORRELATION BY CKB WILL 
+ APPLY, ',/, WITHIN ACCEPTABLE’, 
+' ENGINEERING LIMITS WHEN THESE TEMPERATURES',/, 
+'EXCEED THE MODELS RANGE.'//, "MIXTURE TEMPERATURE IS 
+ ',F7.2,' DEG R',/) 


211 FORMAT (//,6(F10.4,2X)) 
GO TO 260 
Cc 
250 WRITE(7,255) PR,TR,XG1,XG2,T,Z,M 
255 FORMAT (‘ERROR WAS DETECTED’, 6G12.4) 
260 STOP 


END 
ULB Oo Boe 2 Saeco 
@ 
SUBROUTINE MIXT(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10, 
+X11,X12,X13,X14,PC1,TC1,H2S 1,CO21,XN21,PC2,TC2, 
+H2S2,CO22,XN22,PC3,TC3,G1,G2,SM1,SM2) 
DIMENSION TC(14), VC(14), WM(14), ZC(14), ZX(14) 
ZX(1)=X1 
ZX(2)=X2 
ZX (3)=X3 
ZX(4)=X4 
ZX(5)=X5 
ZX(6)=X6 
ZX(7)=X7 
ZX(8)=X8 
ZX(9)=X9 
ZX(10)=X10 
ZX(11)=X11 
ZX(12)=X12 
ZX(13)=X13 
ZX(14)=X14 
C 
DATA TC /343.32, 550.10, 666.01, 765.62, 734.96, 
z 845.60, 829.03, 913.60, 972.70, 1024.1, 
: 547.87, 672.60, 227.30, 9.69/ 
C 
DATA VC /1.5851, 2.3695, 3.2500, 4.0861, 4.2082, 
x 4.8630, 4.8991, 5.9204, 6.9142, 7.8820, 
+ 1.5051, 1.5675, 1.4455, 0.9207/ 
é 
DATA WM /16.043, 30.070, 44.097, 58.124, 58.124, 
= 72.151, 72.151, 86.178, 100.205, 114.232, 
+ 44.010, 34.076, 28.013, 4.003/ 
G 
DATA ZC /0.2876, 0.2844, 0.2805, 0.2742, 0.2826, 
+ 0.2621, 0.2703, 0.2641, 0.2631, 0.2589, 
‘. 0.2745, 0.2252, 0.2925, 0.2923/ 
C 
PTC = 0.0 
SZC = 0.0 
SVC = 0.0 


SM1 = 0.0 
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DO 40 II=1,14 

PTC = PTC + ZX(IJ) * TC(II)/100.0 

SZC = SZC + ZX(ID) * ZC(ID/100.0 

SVC = SVC + ZX(ID * VC(ID/100.0 

SM1 = SM1 + ZX(II) * WM(ID/100.0 
CONTINUE 


CO21 = ZX(11) 
H2S1 = ZX(12) 
XN21 =ZX(13) 


G1=SM1/28.97 
R=10.732 
PPCGSREPTC*SZC/SVC 
Pe2=PLC 

PC2=PPC 


IF(GO21--EQ:;0:0-AND.-H2S1..EQ:-0:0)-GO'TO 51 


CALE WACOR((H2S1,CO21,PPG;7PFE.C,PCITC1) 
ORO 32 


TC1=PTC 
PC1=PPC 

TC3 = 0.0 
SZC = 0.0 
SVC = 0.0 
SM2 = 0.0 


DM=100.0-ZX(11)-ZX(12)-ZX(13) 
ZX(11)=0.0 
ZX(12)=0.0 
ZX(13)=0.0 


DO 53 Il=1,14 

TC3 = TC3 + ZX(I) * TC(I)/DM 

SZC = SZC + ZX(I) * ZC(ID)/DM 

SVC = SVC + ZX(II) * VC(ID)/DM 

SM2 = SM2 + ZX(II) * WM(ID/DM 
CONTINUE 


CO22 = 0.0 
H2S2 = 0.0 
XN22 =0.0 


G2=SM2/28.97 
R=10.732 
PES = R*TC3*SZC/SVC 


RETURN 
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CORRECTION FOR H2S AND CO2 CONTAMINATION 
BY WICHERT AND AZIZ CORRELATION 


SUBROUTINE WACOR (H2S,CO2,PPC,PTC,PC,TC) 


A=(CO2+H2S)/100.0 

B=H2S/100.0 
E3=120.0*(A**0.9-A**1.6)+15.0*(B**0.5-B**4,0) 
TC=PTC-E3 

PC=PPC*TC/(PTC+(B*(1-B)*E3)) 

RETURN 

END 


SUBROUTINE BOUND (TR,TC,G1,G2,ERR,T) 
T=TR*TC 


CHECKING UPPER AND LOWER BOUNDARY AS FOLLOWS 
TR - 1.05 TO 3.00 
677055 TO150 
T - 500.0 TO 860.0 DEG R 


IF OUTSIDE OF THESE RANGES THEN GO TO PROGRAM'S END 


IF (TR .LT. 1.05 .OR. TR .GT. 3.00) ERR=1.0 
IF (G1 .LT. 0.55 .OR. G1 .GT. 1.50) ERR=21.0 
IF (G2 .LT. 0.55 .OR. G2 .GT. 1.50) ERR=22.0 
IF (T .LT. 500.0 .OR. T .GT. 860.0) ERR=3.0 
RETURN 

END 


NATUAL GAS VISCOSITY DETERMINATION USING DEAN AND STIEL 
CORRELATION AS PRESENTED IN AIChE MAY 1965. 


SUBROUTINE VIS1 (ZI, VR,DR,V 1) 


VG = 10.8E-05*(EXP(1.439*DR)-EXP(-1.11*DR**1.858)) 
V1 = 10000.0*(VG+VR)/ZI 

RETURN 

END 
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NATURAL GAS VISCOSITY DETERMINATION USING CARR 
KOBAYASHI AND BURROWS CORRELATION AND LAGRANGE 
INTERPOLATION 


SUBROUTINE VIS2 (PR,TR,G,T,H2S,CO2,XN2,V2) 


DIMENSION TRBL(13), PRBL(22), VBL (22,13), VB1(22,5), 
+ VB1A(22,2), VB2(22,6) 


EQUIVALENCE (VB1(1,1),VBL(1,1)),(VB1A(1,1), VBL(1,6)), 
+ (VB2(1,1), VBL(1,8)) 


DATA TRBL / 

+4 1205. T1910) 1.455.1.2001.30,. 1.40, 1.50, 

+00, 177352 2.00)'2225; 2:50, 3.00/ 

DATA PRBL / 

+1018 0.2" 0.350 490.5, 0.6; 0:7, 0.8, 0.9, 1.0;.1:2 
oregon .G. 18s 2093.00 4:0-6:058.0, 10.0,-15.0; 
+ 20.0/ 


DATA VB1/ 
+1.000,1.012,1.025,1.050,1.075,1.100,1.145,1.195,1.285, 
+1.415,1.760,2.285,2.865,3.290,3.650,4.760,5.500,6.460, 
+7.150,7.680,8.650,9.370, 
+1.000,1.011,1.023,1.043,1.065,1.086,1.120,1.150,1.195, 
+1.255,1.435,1.700,2.070,2.465,2.800,3.850,4.655,5.720, 
+6.500,7.060,8.100,8.880, 
+1.000,1.010,1.021,1.036,1.055,1.073,1.095,1.120,1.145, 
+1.175,1.280,1.420,1.590,1.850,2.160,3.225,3.975,5.030, 
+5.820,6.385,7.410,8.180, 
+1.000,1.009,1.019,1.030,1.045,1.060,1.070,1.085, 1.110, 
+1.135,1.195,1.285,1.425,1.570,1.750,2.600,3.350,4.380, 
+5.125,5.740,6.750,7.500, 
+1.000,1.008,1.017,1.027,1.040,1.054,1.063,1.075,1.100, 
+1.120,1.155,1.215,1.285,1.360, 1.460,2.020,2.560,3.500, 
+4.185,4.755,5.790,6.500/ 


DATA VBI1A/ 
+1.000,1.007,1.015,1.024,1.035,1.048,1.056,1.067,1.089, 
F1t100)10 195,121.85) 1.235)1.280,1.335,1.690,2. 110,2.790, 
+3.380,3. 860,4.790,5.410, 
+1,000,1.006, 1.013, 1.021,1.030,1.042,1.049,1.059, 1.076, 
+1.100,1.120,1.150,1.185,1.220,1.260,1.500,1.785,2.325, 
+2.820,3.230,4.060,4.610/ 


DATA VB2/ 
+1.000,1.005,1.011,1.018,1.025,1.036,1.042,1.051,1.067, 


+1.070,1.095,1.120,1.150,1.180,1.215,1.385,1.595,2.030, 
+2.425,2.770,3.490,4.025, 
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+1.000,1.004,1.009,1.015,1.021,1.030,1.035,1.043,1.056, 
rt -00.5312090, 1 110, 1-125,1:145,1.165,1,280,1.435,1.770, 
+2,095-2.51952.990,3.500, 
+1.000,1 003, 1G 007, i 012, 1.017,1.024,1.028,1.035,1.045, 
+1.055,1.060,1.070,1.080,1.095,1.110,1.205,1.290, 1.500, 
#17 2951.955,2-480'2.925, 
+5. 000, 1.002,1.005,1. 009, 1FO01S,1.018).1.0215 1027-1034) 
+1.040,1.045,1.055,1.065,1.075,1.085,1.145.1.210.1.340, 
+1.485,1.665,2.085,2.460, 
+1.000,1.001,1.003,1.006,1.009,1.012,1.015,1.019,1.023, 
+1.025,1.030,1.040,1.050,1.060,1.065,1.105,1.155,1.245, 
+1.360,1.485,1.830,2.150, 
+1.000,1.000,1.001,1.003,1.005,1.007,1.009,1.011,1.013, 
+1.015,1.020,1.025,1.030,1.035,1.040,1.060,1.085,1.140, 
+19205,1.265,1:495,1.750/ 


INITIALIZING UPPER AND LOWER BOUNDARIES FOR THE 
LOOK-UP TABLES 


LOPR = 2 
LOTR =3 
THIPR = 21 
ITHITR = 11 
IERR = 0 


LOCATE TR IN INDEX 


DO 100 J = LOTR,IHITR 
IF (TRBL(J) .GE. TR) GO TO 120 


100 CONTINUE 


CG 
Cc 
8: 
C 
C 
C 
& 
& 
C 
Cc 
Cc 


J = IHITR+1 
CHECK PR AGAINST INDEX 


IF PR<0.01; PROFORM LAGRANGE INTERPOLATION 
USING RATIO FOR PR=0.01 


IF PR>20.0; PROFORM LAGRANGE INTERPOLATION 
USING RATIO FOR PR=20.0 


OTHERWISE LOCATE PR 


120 IF (PRBL(LOPR-1) .LT. PR) GO TO 160 


I=1 
ICALL = 1 
GO TO 320 


140 VISR = 1.0+(PR*(VISI-1.0)) 


GO TO 300 


160 IF (PRBL(IHIPR+1) .GT. PR) GO TO 200 


I = IHIPR+1 
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TeAWIels 
GO TO 320 
180 VISR = VISI 
GO TO 300 
200 DO 2201 = LOPR,IHIPR 
IF (PRBL(I) .GE. PR) GO TO 240 
220 CONTINUE 
I = IHIPR+1 
240 ICALL =3 
GO TO 320 
260 VISJ = VISI 
T= 
ICALL = 4 
GO TO 320 
280 I=I+1 
VISR = VISI+((PR-PRBL(I-1))/((PRBL(1)-PRBL(I-1)))* 
+  (VISJ-VISI 


¢ 
C CORRECTION FOR H2S, CO2, AND N2 
e 


300 NT = T-460.0 

VISU = 0.126585E-01-0.611823E-02*G+0.164574E-02*G**2.0 
+  +0.164574E-04*NT-0.719221E-06*G*NT-0.609046E-06* 
+  G**2,.0*NT 

CORH2S = (0.113E-03*H2S*G-0.38E-04*H2S+0.10E-05)* 

+ (1.0/(1.0+G))+0.10E-05 

CORCO2 = (0.134E-03*CO2*G-0.40E-05*CO2+0.40E-05*G)* 
+ (1.0/(1.0+G))-0.30E-05 

CORN? = (0.170E-03*KN2*G+0.21E-04*XN2+0.10E-04*G)* 
+ (1.0/(1.0+G))-0.60E-05 

VISA = VISU+CORH2S+CORCO2+CORN2 

V2 = 10000.0*VISR*VISA 


Le 
GO TO 360 

Ge 

C CALL LAGRANGE INTERPOLATION 

G 

320 CALL XLGR4(TR,TRBL(J-2),TRBL(J-1), TRBLUJ),TRBL(VJ+1), 


+  VISI,VBL(,J-2),VBL(LJ-1),VBL(,J), VBL(I,J+1)) 
GO TO (140,180,260,280), ICALL 


SUBROUTINE XLGR4(X,X1,X2,X3,X4,Y,Y1,Y2,Y3,Y4) 


Al = X1-X2 
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A2 = X1-X3 

A3 = X1-X4 

A4 = X2-X3 

AS = X2-X4 

A6 = X3-X4 

Bl =X-X1 

B2 = X-X2 

B3 = X-X3 

B4 = X-X4 

Y = B2/A1*B3/A2*B4/A3*Y 1-B1/A1*B3/A4*B4/A5*Y2+ 
+ B1/A2*B2/A4*B4/A6*Y3-B1/A3*B2/A5*B3/A6* Y4 
RETURN 

END 


C 
c 
C NATURAL GAS VISCOSITY DETERMINATION USING LEE GONZALEZ 
C AND EAKIN CORRELATION AS PRESENTED IN SPE AUG 1966 

Ss 

SUBROUTINE VIS3 (XK,XX,XY,DE,V3) 

C 


V3 = XK*EXP(XX*DE**XY) 
RETURN 


Ce) 


SUBROUTINE ZKATZ(PR,TR,Z,DR,ITER) 


DIMENSION A(8) 
DATA A/0.31506237,-1.0467099,-0.57832729,0.53530771, 
+  -0.6123032,-0.10488813,0.68157001,0.68446549/ 
ITER =0 
DR = 1.0 
IF(TR-1.05) 100, 10,10 
10 IF(TR-3.0)20,20,100 
20 IF(PR-30.0)30,30,100 
30 DO 80 ITER = 1,10 
DR2 = DR**2 
T1 = (A(1)*TR+A(2)+A(3)/TR**2)*DR 
T2 = (A(4)*TR+A(5))*DR2 
T3 = A(5)*A(6)*DR**5 
T4 = A(7)*DR2/TR**2 
T5 = A(8)*DR2 
T6 = EXP(-TS) 
P = (TR+T1+T2+T3)*DR+T4*DR*(1.0+T5)*T6 
DP = TR+2.0*T 14+3.0*T2+6.0*T3+T4*T6*(3.0+3.0* 
+ T5-2.0*T5**2) 
DR1 = DR-(P-0.270*PR)/DP 
IF(DR1)40,40,50 
40 DR1 =0.5*DR 
50 IF(DR1-2.2)70,70,60 
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60 DR1 = DR+0.9*(2.2-DR) 
70 IF(ABS(DR-DR1)-0.1E-04)90,80,80 
80 DR=DR1 
90 Z=0.270*PR/(DR1*TR) 
DR = 0.270*PR/(Z*TR) 
DZ = A(1)+A(2)/TR+A(3)/TR**3 + 2.0*(A(4)+A(5)/TR)*DR + 
+ 5.0*A(5)*A(6)*DR**4/TR + (1.0+A(8)*DR**2- 
+  A(8)**2**DR**4)*2.0*A(7)*DR/TR**3*EXP(-A(8)* 
+ DR**2) 
CR = 1.0/PR -0.270*DZ/(Z**2*TR*(1.0 + DR/Z*DZ)) 
100 RETURN 
END 
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APPENDIX V 


THIS PROGRAM IS USE TO DETERMINE THE PSEUDO GAS 
PRESSURE INTERGRAL USING THE TRAPIZODIAL RULE. 
COMPRESSILITY FACTOR WAS OBTAINED 
USING DRANCHUK, PURVIS, AND ROBERTSON'S 
MODEL OF THE BWR SOLUTION, AND VISCOSITY 
BY DEAN AND STIEL MODEL, AND CARR, KOBAYASHI 
AND BURROWS 


CARR, KOBAYASHI AND BURROWS CORRELATION AND LAGRANGE 
INTERPOLATION WERE PROVIDED BY THE ERCB IN THEIR 
APPENDIX 'D', PROGRAMS 


A COMPARISON OF THE EFFECTS OF CO2,H2S AND N2 WERE ALSO 
MADE BETWEEN THE METHODS WITH THE RESULTS STORED IN 
FILES 6-10. 


A LISTING OF THESE RESULTS HAVE BEEN ATTACHED 
NOMENCLATURE USED WITHIN THIS PROGRAM 


A - MOLAR FRACTION OF H2S + CO2 
B - MOLAR FRACTION OF H2S 
DIFF1-DIFFERENCE IN % BETWEEN OTHER XINTG AND XINTG2 
DIFF2-DIFFERENCE IN % BETWEEN OTHER XINTG AND XINTG3 
DIFF3-DIFFERENCE IN % BETWEEN OTHER XINTG AND XINTG4 
DIFF4-DIFFERENCE IN % BETWEEN OTHER XINTG AND XINTG5 
DIFF5-DIFFERENCE IN % BETWEEN OTHER XINTG AND XINTG6 
E3 - WICHERT & AZIZ CORRECTION FOR ACID GASES 
G -GAS GRAVITY (ENTER 0.0 INDATA READ STATEMENT 
TO CALCULATED G FROM COMPOSITION) 
H_ - STEP SIZE INTERVAL 
XINTG1 - PSEUDO PRESSURE CALCULATED USING CKB'S 
VISCOSITY AND WA CORRELATION TO CORRECT FOR 
ACID GASES 
XINTG2 - PSEUDO PRESSURE CALCULATED USING DS'S 
VISCOSITY AND WA CORRELATION TO CORRECT FOR 
ACID GASES 
XINTG3 - PSEUDO PRESSURE CALCULATED USING CKB'S 
VISCOSITY W/O WA CORRELATION TO CORRECT FOR 
ACID GASES, HOWEVER H2S AND CO2 WERE 
INCLUDED IN DETERMINING CRITICAL PROPERTIES 
XINTG4 - PSEUDO PRESSURE CALCULATED USING CKB'S 
VISCOSITY W/O ANY H2S OR CO2 EFFECTS. 
XINTGS5 - PSEUDO PRESSURE CALCULATED USING DS'S 
VISCOSITY W/O WA CORRELATION TO CORRECT FOR 
ACID GASES, HOWEVER H2S AND CO2 WERE 
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INCLUDED IN DETERMINING CRITICAL PROPERTIES 
XINTG6 - PPEUDO PRESSURE CALCULATED USING DS's 
VISCOSITY W/O ANY H2S OR CO2 EFFECTS. 

N - NUMBER OF INTERVAL STEPS 

NN - NUMBER OF PSEUDO-PRESSURE TO CALCULATE 

NP - CRITICAL PRESSURE ATM 

NT - CRITICAL TEMPERATURE K 

PPC - CRITICAL PRERRURE PSIA 

PR - REDUCED PRESSURE 

PX1:4- PRESSURE STEP IN THE TRAPIZODIAL INTEGRAL 

Q - FLOW RATES MMSCF/D 

R_ - UNIVERSAL GAS CONSTANT 
10.732 PSIA-CU-FT/DEG-R-LBMOL 

S1 - SWITCH-ELIMINATING CO2,H2S&N2 FROM GAS 
COMPOSITION 

S2 - SWITCH-ELIMINATES CO2,H2S,N2 FROM CKB MODEL 

S3 - SWITCH-ELIMATES C02&H2S FROM WA CORRECTION 

SM - PSEUDO MOLECULAR WEIGHT LB/LBMOL 

T -SANDFACE TEMPERATURE DEG R 

TC - CRITICAL TEMPERATURE DEGR 

TCO - OLD CRITICAL TEMPERATURE FOR USE IN WICHERT & 
AZIZ CORRELATION 

TI - TIME hrs 

TP - PREVIOUS PRODUCING TIME 

TR - REDUCED TEMPERATURE 

TT - CRITICAL TEMP DEG R 

VC - CRITICAL VOLUME CU-FT/LBMOL 

V -GAS VISCOSITY AS CALCULATED BY DEAN & STIEL 
METHOD; cp 

VIS2- GAS VISCOSITY AS CALCULATED BY CARR, KOBAYASHI 
& BURROWS; cp 

VISR- VISCOSITY RATIO FROM CARR, KOBAYASHI & BURROWS 
MODEL 

VG - VISCOSITY FACTOR LESS THE VISCOSITY EFFECTS AT 
MODERATE PRESSURE (0.2-5 ATM) 

VR - RESIDUAL VISCOSITY 

WM - MOLECULAR WEIGHT LB/LBMOL 

XP - CRITICAL PRESSURE ATM 

XT - HORNER'S TIME RATIO 

ZC - CRITICAL COMPRESSIBILITY FACTOR 

ZI - VISCOSITY PARAMETER 


PROPERTIES LISTED BELOW HAVE BEEN OBTAINED FROM 
GAS PROCESSOR SUPPLIERS ASSOCIATION, ENGINEERING 
DATA BOOK, NINTH EDITION 


CONVERSION USED TO OBTAINED UNITS USED IN THE PROGRAM 
ARE: 


TC = TC(OF) + 460.0 
VC = VC(FT3/LB)*WM(LB/LBMOL) 
ZC = PC*VC/R*TC 
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GAS COMPOSITION IF KNOWN, SHOULD BE EXPRESSED IN 
PERCENTAGES AND ENTERED IN THE FOLLOWING ORDER 


Oly C2NG3 nC) 1C4, n€5,1€5,C6-C7 
C+, 'CO2; H2S, N2, He 


CRITICAL PROPERTIES HAVE BEEN DETERMINED BY KAY'S RULE 


EXCEPT FOR CRITICAL PRESSURE WHICH HAS BEEN 
DETERMINED BY PRAUSNITZ AND GUNN MIXING RULE. 


(oT OLGLG (Pe L@OLOL eS (1 CLOG Sie 


DIMENSION ZPC(5),ZTC(5),ZERR(5),ZF(5),ZT (5), TT (5) 
DIMENSION XP(5),X(14), TP(5),ZP(6),F(6) 
DOUBLE PRECISION ZX1,ZX2,ZX3,VC,SM,SVC,TT,XP,V1,V2,Z, 
DOUBLE PRECISION G1,G2 
TR=1.05 
READ (5,5) NN 
5 FORMAT (12) 
WRITE (6,10) TR 
WRITE (7,11) TR 
10 FORMAT (2X,'THE FOLLOWING CALCULATIONS ARE FOR A', 
+'REDUCED TEMPERATURE OF |,F4.2,//, 
+ 2X,'PSEUDO-PRESSURE USING VISCOSITY DETERMINED’ 
+,/,2X,! BY CARR KOBAYASHI BURROWS',/, 
+/,'MIXTURE',3X,'PR',5X,'PRESSURE’, 
+4X,'PSIA2/up',2X,'%DIFF',/, 19X, 
+ 'PSIA',/) 
FORMAT (2X,'THE FOLLOWING CALCULATIONS ARE FOR A’, 
+'REDUCED TEMPERATURE OF |,F4.2,//, 
+ 2X,'PSEUDO-PRESSURE USING VISCOSITY DETERMINED' 
+,/,2X,! BY DEAN & STIEL'/, 
+/,'MIXTURE',3X,'PR',3X,' PRESSURE’, 
+4X,'PSIA2/up',2X,'%DIFF',2X,'%DIFF',2X,'%DIFF',/,19X, 
WUPSTARO4AXI 1-2"3X11/53X,1-6') 


1 


p= 


DO 240 I=1,NN 
20 READ (5,30) (X(J), U=1,14) 
30 FORMAT (14(F5.2,1X)) 


CALL MIXT(X(1),X(2),X(3),X(4),X(5),X(6),X(7),X(8) 
+,X(9),X(10),X(11),X(12),X(13),X(14),ZPC(1),ZTC(1), 
+H251,C021,N21,ZPC(2),ZTC(2),H2S2,C022,N22, 
+ZPC(3),ZTC(3),XG1,XG2,SM1,SM2) 


WRITE (8,32) I,ZPC(1),ZTC(1),ZPC(2),ZTC(2),ZPC(3) 
+ ,ZTC(3),SM1,XG1,SM2,XG2 

32. FORMAT (12,2X,12(F9.3,2X)) 
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CHECKING UPPER AND LOWER BOUNDARY AS FOLLOWS 


DO 35 IA=1,3 

ZERR(IA)=0.0 

CALL BOUND (TR,ZTC(IA),XG1,XG2,ZERR(IA),ZT(IA)) 
35 CONTINUE 
C 


CTS) 


PR=0.25 


C NATUAL GAS VISCOSITY FACTORS FOR THE DETERMINATION OF 
C DEAN AND STIEL VISCOSITY CORRELATION 
IF(TR-1.5) 80,90,90 
80 VR=34.0E-05*TR**(8.0/9.0) 
GO TO 100 
90 VR=166.8E-05*(0.1338*TR-0.0932)**(5.0/9.0) 
100 DO 105 IC=1,3 
TT(IC)=ZTC(IC)*5.0/9.0 
XP(IC)=ZPC(IC)/14.7 
105 CONTINUE 
ZX 1=TT (1)**(1.0/6.0)/(SM1**0.5*XP(1)**(2.0/3.0)) 
ZX2=TT (2)**(1.0/6.0)/(SM1**0.5*XP(2)**(2.0/3.0)) 
ZX3=TT(3)**(1.0/6.0)/(SM2**0.5*XP(3)**(2.0/3.0)) 


DO 205 IB=1,50 
IF (PR .LT. 0.01) GO TO 160 
IF (PR .GT. 20.0) GO TO 220 


CALL ITGR1 (ZPC(1),ZTC(1),PR,TR,XG1,ZT(1), 
+H2S1,CO21,N21,—FLAG1,ZPC(1),XINTG1) 
CALL ITGR2 (PR,TR,VR,ZX1,ZPC(1),XINTG2) 
CALL ITGR1 (ZPC(2),ZTC(2),PR,TR,XG1,ZT (2), 
+H2S1,CO21,N21,FLAG2,ZPC(2),XINTG3) 
CALL ITGR1 (ZPC(3),ZTC(3),PR,TR,XG2,ZT (3), 
+H2S2,CO22,N22,FLAG3,ZPC(3),XINTG4) 
CALL ITGR2 (PR,TR,VR,ZX2,ZPC(2),XINTGS5) 
CALL ITGR2 (PR,TR,VR,ZX3,ZPC(3),XINTG6) 


DIFF1 = (XINTG1-XINTG2)*100.0/XINTG1 
DIFF2 = (XINTG1-XINTG3)*100.0/XINTG1 
DIFF2A= (XINTG2-XINTG3)*100.0/XINTG2 
DIFF3 = (XINTG1-XINTG4)*100.0/XINTG1 
DIFF3A= (XINTG2-XINTG4)*100.0/XINTG2 
DIFF3B= (XINTG3-XINTG4)*100.0/XINTG3 
DIFF4 = (XINTG1-XINTG5)*100.0/XINTG1 
DIFF4A= (XINTG2-XINTG5)*100.0/XINTG2 
DIFF4B= (XINTG3-XINTGS5)*100.0/XINTG3 
DIFF4C= (XINTG4-XINTGS5)*100.0/XINTG4 
DIFF5 = (XINTG1-XINTG6)*100.0/XINTG1 
DIFF5A= (XINTG2-XINTG6)*100.0/XINTG2 
DIFF5B= (XINTG3-XINTG6)*100.0/XINTG3 
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DIFF5C= (XINTG4-XINTG6)*100.0/KINTG4 
DIFFSD= (XINTGS-XINTG6)*100.0/XINTG5 


DO 140 IC=1,3 
ZP(IC) = PR*ZPC(IC) 
140 CONTINUE 
& 


G 


WRITE(6,155) I,PR,ZP(1), XINTG1,XINTG2,XINTG3,XINTG4 
+,XINTG5,XINTG6,DIFF1,DIFF2,DIFF3,DIFF4,DIFF5 
‘G 
WRITE(7,155) I,PR,ZP(1), DIFF2A,DIFF3A,DIFF4A, 
+DIFF5A,DIFF3B,DIFF4B,DIFF5B,DIFF4C,DIFF5C,DIFF5D 
155 FORMAT(2X,12,4X,F5.2,3X,F8.1,2X,11(F10.4,3X)) 
ee 


IF (PR-0.29) 160,165,165 
160 PR=PR+0.05 

GO TO 205 
165 IF (PR-1.49) 170,175,175 
170 PR=PR+0.10 


1752.1Fi(PR-2.99) 180,185,185 
180 PR=PR+0.50 

H=0.1 

GO TO 205 
1857 1Fi(PR-9.99))190,195,195 
190 PR=PR+1.0 


195 IF(PR-19.999) 200,240,240 
200 PR=PR+5.0 
H=0.5 
205 CONTINUE 
220 WRITE (6,212) 
WRITE (7,212) 
212 FORMAT (''/) 
240 CONTINUE 


WRITE (6,211) ZTC(1),ZTC(2),ZTC(3),ZPC(1),ZPC(2), 
+ZPC(3) 
WRITE (7,210) ZT(1) 
210 FORMAT (//,,TEMPERATURE IS OUTSIDE RANGE OF CKB', 
+! MODEL WHERE TEMP>860:'/,,1T IS ASSUMED THAT 
+ VISCOSITY',' CORRELATION BY CKB WILL 
+ APPLY, ',/,,;WITHIN ACCEPTABLE’, 
4' ENGINEERING LIMITS WHEN THESE TEMPERATURES',/, 
4'EXCEED THE MODELS RANGE.’ ‘MIXTURE TEMPERATURE IS 
+ '\F7.2,' DEG R',/) 
211 FORMAT (//,6(F10.4,2X)) 
GO TO 260 
e 
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250 WRITE(7,255) PR,TR,XG1,XG2,T,Z,M 
255 FORMAT (ERROR WAS DETECT ED', 6G12.4) 
260 STOP 

END 


SUBROUTINE MIXT(X1,X2,X3,X4,X5,X6,X7,X8,X9,X10, 
+X11,X12,X13,X14,PC1,TC1,H2S1,CO21,N21,PC2,TC2, 
+H2S2, CO22, N22 G31 E2.GL G2. SM1 ,SM2) 


DIMENSION TC(14), VC(14), WM(14), ZC(14), ZX(14) 


REAL N2 
ZX(1)=X1 
ZX (2)=X2 
ZX(3)=X3 
ZX(4)=X4 
ZX(5)=X5 
ZX(6)=X6 
ZX(7)=X7 
ZX(8)=X8 
ZX(9)=X9 
ZX(10)=X10 
ZX(11)=X11 
ZX(12)=X12 
ZX(13)=X13 
ZX(14)=X14 


DATA TC /343.32, 550.10, 666.01, 765.62, 734.96, 
+ 845.60, 829.03, 913.60, 972.70, 1024.1, 
+ 5478 75072.60, 227.30; 9.69/ 


DATA VC /1.5851, 2.3695, 3.2500, 4.0861, 4.2082, 
+ 4.8630, 4.8991, 5.9204, 6.9142, 7.8820, 
+ 1.5051, 1.5675, 1.4455, 0.9207/ 


DATA WM /16.043, 30.070, 44.097, 58.124, 58.124, 
+ T2ava), 72151, 86.178, 100.205; 114.232, 
+ 44.010, 34.076, 28.013, 4.003/ 


DATA ZC /0.2876, 0.2844, 0.2805, 0.2742, 0.2826, 
+ 0.2621, 0.2703, 0.2641, 0.2631, 0.2589, 
+ 027457 0:2252;,0.2925, 0.2923/ 


PTC = 0.0 
SZC = 0.0 
SVC = 0.0 
SM1 = 0.0 


DO 40 II=1,14 
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PTC = PTC + ZX(II) * TC(I1)/100.0 

SZC = SZC + ZX(II) * ZC(ID/100.0 

SVC = SVC + ZX(II) * VC(ID/100.0 

SM1 = SMI + ZX(I) * WM(ID/100.0 
CONTINUE 


CO21 = ZX(11) 
H2S1 = ZX(12) 
N21 =ZX(13) 


G1=SM1/28.97 
R=10.732 

PPC = R*PTC*SZC/SVC 
TC2=PTC 

PC2=PPC 


TCO? Tt EOs0.0.AND, 2S 1. EO:0:0).GO TO51 


GALL WACOR (H2S81,C021,PPC,PTC,PC1,TC1) 
GO TO 52 


TC1=PTC 
PC1=PPC 

TC3 = 0.0 
SZC = 0.0 
SVC = 0.0 
SM2 = 0.0 


DM=100.0-ZX(11)-ZX(12)-ZX(13) 
ZX(11)=0.0 
ZX(12)=0.0 
ZX(13)=0.0 


DO 53 II=1,14 

TC3 = TC3 + ZX(I) * TC(ID)/DM 

SZC = SZC + ZX(Il) * ZC(I)/DM 

SVC = SVC + ZX(II) * VC(ID/DM 

SM2 = SM2 + ZX(ID) * WM(ID/DM 
CONTINUE 


CO22 = 0.0 
H2S2 = 0.0 
N22 =0.0 


G2=SM2/28.97 
R=10.732 
PC3= R*T.C3*SZC/SVC 


RETURN 
END 
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CORRECTION FOR H2S AND CO2 CONTAMINATION 
BY WICHERT AND AZIZ CORRELATION 


SUBROUTINE WACOR (H2S,CO2,PPC,PTC,PC,TC) 
A=(CO2+H28)/100.0 


B=H2S/100.0 
E3=120.0*(A**0.9-A**1.6)+15.0*(B**0.5-B**4.0) 
TC=PTC-E3 
PC=PPC*TC/(PTC+(B*(1-B)*E3)) 
RETURN 
END 
SUBROUTINE BOUND (TR,TC,G1,G2,ERR,T) 
T=TR*TC 
CHECKING UPPER AND LOWER BOUNDARY AS FOLLOWS 
TR - 1.05 TO 3.00 
Gr 0)95: 107250 
T - 500.0 TO 860.0 DEG R 
IF OUTSIDE OF THESE RANGES THEN GO TO PROGRAM'S END 
IF (TR .LT. 1.05 .OR. TR .GT. 3.00) ERR=1.0 
IF (G1 .LT. 0.55 .OR. G1 .GT. 1.50) ERR=21.0 
1H (G2e6770.55 .OR. G2 .GT. 1.50) ERR=22.0 
IF (T .LT. 500.0 .OR. T .GT. 860.0) ERR=3.0 
RETURN 
END 


SUBROUTINE ITGR1 (PC1,TC1,PR,TR,G1,T1, 
+H2S1,CO21,N21,FLAG1,PC,XINTG1) 


SUM1=0.0 
PR1=0.20 
H=0.05 
PX1A=0.0 


THE FOLLOWING INTERGRALS ARE DONE BY 
TRAPIZODIAL 


DO 145 II=1,400 
IF (PR1 .GT. PR) GO TO 150 
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CALL ZKATZ (PR1,TR,Z,DR,ITER) 
IF (ITER .EQ. 0) FLAG1=0.0 
CALL VIS2 (PR1,TR,G1,T1,H2S1,CO21,N21,V2) 
PX1 = 2.0*(PR1/(Z*V2)) 
SUM1 = SUM1 + H*(PX1 + PX1A)/2.0 
PX1A = PX1 
PR1=PR1+H 

145 CONTINUE 

C 

150 XINTG1 = SUM1*PC**2.0/1000.0 
RETURN 
END 


SUBROUTINE ITGR2 (PR,TR,VR,ZZI,PC2,XINTG2) 


H=0.05 
SUM2 = 0.0 
PX2A = 0.0 
PR2 = 0.2 


THE FOLLOWING INTERGRALS ARE DONE BY 
TRAPIZODIAL 


QAO 


DO 145 II=1,400 
IF (PR2.GT.PR) GO TO 150 
CALL ZKATZ (PR2,TR,Z,DR,ITER) 
CALL VIS1 (ZZI,VR,DR,V1) 
PX2 = 2.0*(PR2/(Z*V1)) 
SUM2 = SUM2 + H*(PX2 + PX2A)/2.0 
PX2A = PX2 
PR2=PR2+H 
145 CONTINUE 
C 


150 XINTG2 = SUM2*PC2**2.0/1000.0 
RETURN 
END 


€ 
= 
C NATUALGAS VISCOSITY DETERMINATION USING DEAN AND STIEL 
C CORRELATION AS PRESENTED IN AIChE MAY 1965. 
Gs 
Cc 


SUBROUTINE VIS1 (ZI, VR,DR,V1) 


VG = 10.8E-05*(EXP(1.439*DR)-EXP(-1.11*DR**1.858)) 
V1 = 10000.0*(VG+VR)/ZI 

RETURN 

END 
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NATURAL GAS VISCOSITY DETERMINATION USING CARR 
KOBAYASHI AND BURROWS CORRELATION AND LAGRANGE 
INTERPOLATION 


SUBROUTINE VIS2 (PR,TR,G,T,H2S,CO2,N2,V2) 


DIMENSION TRBL(13), PRBL(22), VBL(22,13), VB1(22,5), 
+ VB1A(22,2), VB2(22,6) 


EQUIVALENCE (VB1(1,1),VBL(,1)),(VB1A(1,1), VBL(1,6)), 
+ (VB2(1,1), VBL(1,8)) 


DATA TRBL/ 

ele yall OF, 1915,01,2031.30, 1.40, 1.50, 
te lnO0 bo.e2 004 225,.2100,3.00/ 
DATA PRBL / 

TaUsi mow, 0.570.4e0.5, 0.6, 0.7,.0.8, 0.9, 1.0; 1.2 
ee Ae Oe bis, 910,74,05.0,0,.8-05.10,0,519.0; 
+ 20.0/ 


DATA VB1/ 
+1.000,1.012,1.025,1.050,1.075,1.100,1.145,1.195,1.285, 
+1.415,1.760,2.285,2.865,3.290,3.650,4.760,5.500,6.460, 
+7.150,7.680,8.650,9.370, 
+1.000,1.011,1.023,1.043,1.065,1.086,1.120,1.150,1.195, 
+1.255,1.435,1.700,2.070,2.465,2.800,3.850,4.655,5.720, 
+6.500,7.060,8.100,8.880, 
+1.000,1.010,1.021,1.036,1.055,1.073,1.095,1.120,1.145, 
+ h17 5512280, 1.420,1.590,1:850,2:160,3.225,3.975,5.030, 
+5.820,6.385,7.410,8.180, 
+1.000,1.009, 1.019, 1.030,1.045,1.060,1.070,1.085,1.110, 
+1.135,1.195,1.285, 1.425,1.570,1.750,2.600,3.350,4.380, 
+5.125,5.740,6.750,7.500, 
+1.000,1.008,1.017,1.027,1.040,1.054,1.063,1.075,1.100, 
+1.120,1.155,1.215,1.285,1.360,1.460,2.020,2.560,3.500, 

+4.185,4.755,5.790,6.500/ 


DATA VBIA/ 
+1.000,1.007,1.015,1.024,1.035,1.048,1.056,1.067,1.089, 
+1.100,1.135,1.185, 1.235,1.280,1.335,1.690,2.110,2.790, 

+3.380,3.860,4.790,5.410, 
+1.000,1.006, 1.013, 1.021,1.030,1.042,1.049,1.059, 1.076, 
+1.100,1.120,1.150, 1.185, 1.220, 1.260, 1.500, 1.785,2.325, 
+2.820,3.230,4.060,4.610/ 


DATA VB2/ 
+1.000,1.005,1.011,1.018,1.025,1.036,1.042,1.051,1.067, 
+1 070,1.095, 1.120, 1.150,1.180,1.215,1.385,1.595,2.030, 
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+2.425,2.770,3.490,4.025, 
+1.000,1.004,1.009,1.015,1.021,1.030,1.035,1.043, 1.056, 
ty09316090;1.110,1.125,17145,1.165,1.280,1.435,1.770, 
#21095;2.37532.990,3.500, 
t°000F12003 1.007; 1.012,1.017,1.024,1.028,1.035,1.045, 
+1.055,1.060,1.070, 1.080,1.095,1.110,1.205, 1.290, 1.500, 
af 72931 955,2548032.925; 
+1.000,1.002,1.005,1.009,1.013,1.018,1.021,1.027, 1.034, 
+1.040,1.045,1.055,1.065,1.075,1.085,1.145,1.210,1.340, 
+1.485,1.665,2.085,2.460, 
+1.000,1.001,1.003, 1.006,1.009,1.012,1.015,1.019, 1.023, 
+1.025,1.030,1.040,1.050,1.060,1.065,1.105,1.155,1.245, 
+1.360,1.485,1.830,2.150, 
+1.000,1.000,1.001,1.003,1.005,1.007,1.009,1.011,1.013, 
+1.015,1.020,1.025,1.030,1.035,1.040,1.060,1.085,1.140, 
#1205. .265,1.495,1.750/ 


INITIALIZING UPPER AND LOWER BOUNDARIES FOR THE 
LOOK-UP TABLES 


LOPR = 2 
LOTR = 3 
ITHIPR = 21 
IHITR = 11 
IERR = 0 


LOCATE TR IN INDEX 


DO 100 J = LOTR, IHITR 
IF (TRBL(J) .GE. TR) GO TO 120 


100 CONTINUE 


Cc 
C 
c 
iS 
C 
¢ 
S 
( 
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C 
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J = IHITR+1 
CHECK PR AGAINST INDEX 


IF PR<0.01; PROFORM LAGRANGE INTERPOLATION 
USING RATIO FOR PR=0.01 


IF PR>20.0; PROFORM LAGRANGE INTERPOLATION 
USING RATIO FOR PR=20.0 


OTHERWISE LOCATE PR 


120 IF (PRBL(LOPR-1) .LT. PR) GO TO 160 


ICALL = 1 
GO TO 320 


140 VISR = 1.0+(PR*(VISI-1.0)) 


GO TO 300 


160 IF (PRBL(IHIPR+1) .GT. PR) GO TO 200 
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I = IHIPR+1 
IGAETX=2 
GO TO 320 
180 VISR = VISI 
GO TO 300 
200 DO 2201=LOPR,IHIPR 
IF (PRBL(I) .GE. PR) GO TO 240 
220 CONTINUE 
I = IHIPR+1 
240 FICARLI=3 
GO TO 320 
260 VISJ = VISI 
Fail 
ICALRL =4 
GO TO 320 
280 I=I+1 
VISR = VISI+((PR-PRBL(I-1))/(PRBL(I)-PRBL(I-1)))* 
+ (VISJ-VISD 


Cc 
C CORRECTION FOR H2S, CO2, AND N2 
C 


300 NT = T-460.0 

VISU = 0.126585E-01-0.611823E-02*G+0.164574E-02*G**2.0 
+  +0.164574E-04*NT-0.719221E-06*G*NT-0.609046E-06* 
+  G**2.0*NT 

CORH2S = (0.113E-03*H2S*G-0.38E-04*H2S+0.10E-05)* 

+  (1.0/(1.0+G))+0.10E-05 

CORCO2 = (0.134E-03*CO2*G-0.40E-05*CO2+0.40E-05*G)* 
+ (1.0/(1.0+G))-0.30E-05 

CORN? = (0.170E-03*N2*G+0.21E-04*N2+0.10E-04*G)* 
+ (1.0/(1.0+G))-0.60E-05 

VISA = VISU+CORH2S+CORCO2+CORN2 

V2 = 10000.0*VISR*VISA 


(@; 
GO TO 360 

C 

C CALL LAGRANGE INTERPOLATION 

c 

32 


0 CALL XLGR4(TR,TRBL(J-2), TRBLGU-1),TRBLVJ),TRBLU+1), 


+ VISI,VBL(,jJ-2), VBL(,J-1), VBL(,J), VBL(I,J+1)) 
GO TO (140,180,260,280), ICALL 
340 IERR=1 


SUBROUTINE XLGR4(X,X1,X2,X3,X4,Y,Y1,Y2,Y3,Y4) 
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(int) LUE ees 


G4 


Al = X1-X2 

A2 = X1-X3 

A3 = X1-X4 

A4 = X2-X3 

AS = X2-X4 

A6 = X3-X4 

Bl =xX-X1 

B2 = X-X2 

B3 = X-X3 

B4 = X-X4 

Y = B2/A1*B3/A2*B4/A3*Y 1-B1/A1*B3/A4*B4/A5*Y2+ 
+ B1/A2*B2/A4*B4/A6*Y3-B 1/A3*B2/A5*B3/A6* Y4 
RETURN 

END 


SUBROUTINE ZKATZ(PR,TR,Z,DR,ITER) 


DIMENSION A(8) 
DATA A/0.31506237,-1.0467099,-0.57832729,0.53530771, 
+ -0.6123032,-0.10488813,0.68 157001,0.68446549/ 
ITER=0 
DRO 
IF(TR-1.05)100,10,10 
10 IF(TR-3.0)20,20,100 
20 IF(PR-30.0)30,30,100 
30 DO 80 ITER = 1,10 
DRZ = DR**2 
T1 = (A(1)*TR+A(2)+A(3)/TR**2)*DR 
T2 = (A(4)*TR+A(5))*DR2 
T3 = A(5)*A(6)*DR**5 
4 = A(7)*DR2/TR**2 
T5 = A(8)*DR2 
T6 = EXP(-TS) 
P = (TR+T1+T2+T3)*DR+T4*DR*(1.0+T5)*T6 
DP = TR+2.0*T 1+3.0*T2+6.0*T3+T4*T6*(3.0+3.0* 
‘foe DP o=2.07 bon 2) 
DR1 = DR-(P-0.270*PR)/DP 
IF(DR1)40,40,50 
40 DR1=0.5*DR 
50 IF(DR1-2.2)70,70,60 
60 DR1 = DR+0.9*(2.2-DR) 
70 IF(ABS(DR-DR1)-0.1E-04)90,80,80 
80 DR=DRI1 
90 Z =0.270*PR/(DR1*TR) 
DRe= 0.270*PR/(Z* TR) 
DZ = A(1)+A(2)/TR+A(3)/TR**3 + 2.0*(A(4)+A(5)/TR)*DR + 
+ 5,0*A(5)*A(6)*DR**4/TR + (1.0+A(8)*DR**2- 
+ A(8)**2**DR**4)*2.0*A(7)*DR/TR**3*EXP(-A(8)* 
+ DR**?2) 
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@R=110/PR-0:270*DZ/(Z5*2* TR*(1.0 + DR/Z*DZ)) 
100 RETURN 
END 
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